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ABSTRACT 

Black hole (BH) accretion flows and jets are qualitatively affected by the presence of 
ordered magnetic fields. We study fully three-dimensional global general relativistic magne- 
tohydrodynamic (MHD) simulations of radially extended and thick (height H to cylindrical 
radius R ratio of \H/R\ ~ 0.2-1) accretion flows around BHs with various dimensionless spins 
(a/M, with BH mass M) and with initially toroidally-dominated (0-directed) and poloidally- 
dominated (R—z directed) magnetic fields. Firstly, for toroidal field models and BHs with high 
enough \a/M\, coherent large-scale (i.e. » H) dipolar poloidal magnetic flux patches emerge, 
thread the BH, and generate transient relativistic jets. Secondly, for poloidal field models, 
poloidal magnetic flux readily accretes through the disk from large radii and builds-up to a nat- 
ural saturation point near the BH. While models with \H/R\ ~ 1 and \a/M\ < 0.5 do not launch 
jets due to quenching by mass infall, for sufficiently high \a/M\ or low \H/R\ the polar mag- 
netic field compresses the inflow into a geometrically thin highly non-axisymmetric "magnet- 
ically choked accretion flow" (MCAF) within which the standard linear magneto-rotational 
instability is suppressed. The condition of a highly-magnetized state over most of the horizon 
is optimal for the Blandford-Znajek mechanism that generates persistent relativistic jets with 
> 100% efficiency for \a/M\ > 0.9. A magnetic Rayleigh-Taylor and Kelvin-Helmholtz un- 
stable magnetospheric interface forms between the compressed inflow and bulging jet magne- 
tosphere, which drives a new jet-disk quasi-periodic oscillation (JD-QPO) mechanism. The 
high-frequency QPO has spherical harmonic \m\ = 1 mode period of t ~ 70GM/c 3 for 
a/M ~ 0.9 with coherence quality factors Q > 10. Overall, our models are qualitatively 
distinct from most prior MHD simulations (typically, \H/R\ <K 1 and poloidal flux is lim- 
ited by initial conditions), so they should prove useful for testing accretion-jet theories and 
measuring a/M in systems such as SgrA* and M87. 

Key words: accretion, accretion discs, black hole physics, hydrodynamics, (magnetohydro- 
dynamics) MHD, methods: numerical, gravitation 



1 INTRODUCTION 

Modern black hole (BH) accretion disk theory suggests that angu- 
lar momentum transport, in the past modelled using an o--viscosity 
l |Shakura & Sunyaev|1973"]|Novikov & Thorne|1973[|Thorne|1974| >, 
is due to magnetohydrodynamic (MHD) turbulence driven by the 
magneto-rotational instability (MRI) within a differentially rotating 
disk ( |Balbus &"Ti awley 1991 19981 and due to large-scale mag- 
netic torques within the plunging region of a BH (Znajek||l976| 



MacDonald 1984; Gammie 1999, 


Krolik 1999( |Agol & Krolik 


2000 ; Li 2000 ; Li & Paczynski 2000 , Li 2002 , Reynolds et al. 2006 


Oda et al. 2009 2010, Penna et al. 


2012|l. However, even outside 
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the plunging region, transport may occur via large-scale magnetic 
field threading the disk (e.g. , |Blandford| 1 976[ [Lovelace| 1 976| |Caol 
|201 1[ > or even by ordered field within the disk far beyond the plung- 
ing region (Bisnovatyi-Kogan & Ruzmaiki n| 1 974| |1976| |Narayan| 
|et al.|2003||Bisnovatyi-Kogan & Lovelace|2007^ . 

How are strong, large-scale, or ordered magnetic fields gen- 
erated in BH accretion flows? The MHD turbulence increases the 
magnetic field's strength up to some non-linear saturation point, 
while the solenoidal constraint in electrodynamics that V ■ B = 
for magnetic field B ensures that, within some radius ro, the mag- 
netic field's flux (e.g. J[ " B z dA at the equator for vertical field B : 
and equatorial area A) remains fixed unless added at r . 

Patches of constant polarity magnetic flux might form stochas- 
tically out of an MHD dynamo (Tout & Pringle 1996) as occurs for 
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the toroidal (^-directed) magnetic field as seen in MHD simulations 
l |Davis et al. 2010; Guan & Gammie 2011) and for the large-scale 
poloidal (Tf-z-directed) magnetic field as seen in the Sun ( Babcock| 
|1959| >. However, so far no evidence exists for the development of 
large-scale poloidal field from only a dominant toroidal field ( |De| 
|Villiers et"aT[|2005| >. The MRI-driven MHD turbulence-generated 
poloidal field remains small-scale on the order of the disk height 
(often used to estimate jet field strength, e.g. in Meier 2001 and 
|Livio et al.|2 003 1 and does not sustain relativistic jets due to persis- 
tent mass-loading (Beckwith et al. 2008a). Patches of constant po- 
larity flux may also develop by the Poynting-Robertson drag (PRD) 
effect ( |Contopoulos & Kazanas[|1998), but the PRD effect might 
saturate at low field strengths (Bisnovat yi-Kogan et al.|2002] >. 

An alternative is that ordered constant polarity magnetic flux 
comes from large radii, as implicit in outflow models (Bland- 
ford||1976[|Brandford & Znajekp977l |Blandford & Payne||1982f 



Narayan et al.|2007|[Tchekhovskoy et a l. 2009 2010a b), but trans- 
port of flux from large radii may be inefficient fLubow et al .]1994| 
Spruit & Uzdensky||2005j |Rothstein & Lovelace||2008| |Beckwifh| 
et al.|2009j >. Such magnetic fields can lead to angular momentum 



transport within the disk and out of a rotating BH (the latter having 
some observational basis; Narayan & McClintock 2012). Indeed, 
production of persistent relativistic jets may require steady large- 
scale dipolar (i.e. sub-dominant higher multipoles) fields near BHs 
( |Beckwith et al.|2008a| [McKinney & Blandford|20"09"[ >. 

Observations do show patches of coherent magnetic flux (O) 
surrounding astrophysical system s that can feed BHs. The inter 
stellar medium has O ~ 0.1pc 2 G ( |Lang et al.|l999[|Vallee|201 1 



threads near the Galactic Nucleus have fl) > 0.01pc 2 G {LaRosa 



|et al.|2004[rFerriere|2009]|Vallee|201 j) , and each O-type star with 
a dipolar field with strength B ~ 1000G feeding SgrA* provides up 
to <f> ~ 10~'°pc 2 G. Such a coherent flux might be available for ac- 
cretion near SgrA*, M87, and other active galactic nuclei (AGNs). 
Indeed, the constancy over several years of the sign of the circular 
polarization from near the BH in SgrA* implies that the magnetic 
field is coherent with constant polarity for many dynamical times 
l |Munoz et al.|20i2) . Also, many AGN jets show persistent signs 
for the transverse Faraday rotation gradient, indicating accretion of 
a persistent magnetic polarity (e.g., Mahmud et al. 2009; Broder- 
|ick & McKinney 2010). For BH x-ray binaries, some portion of 
the donor star's surface dipolar field of order B ~ 100-1000G 
may be available for accretion, such that O ~ 10~ 13 -10~ 12 pc 2 G 
(Justham et al. 2006). For cosmological gamma-ray bursts (GRBs) 
that require B ~ 3 x 10 I5 G for MHD-driven jets by the Blandford- 
Znajek (BZ) mechanism (Blandford & Znajek 1977 , Narayan et al. 
1 1992} |Barkov &~ Komissaro v|2008| , one obtains a requirement of 
<S> ~ 10 9 pc 2 G corresponding to about 10% of the total ordered 
poloidal flux within a presupernova progenitor ( Heger et al .|2005| 
|Komissarov & B arkov 2009 ). Even if A' random polarity patches 
accrete over time, on average 5> oc N ' 2 in a random walk near the 
BH. Some flux can annihilate before reaching the BH, but field re- 
versals are unlikely to have exactly equal flux magnitude and so 
magnetic flux should tend to accumulate. 

Magnetic flux is conserved such that accumulation of a suffi- 
cient amount of constant polarity poloidal magnetic flux brought in 
from large radii might substantially modify the standard picture of 
an MRI-driven MHD turbulent disk that applies for relatively weak 
magnetic field. Accumulation of magnetic flux can eventually lead 
to the formation of a semi-permeable magnetic barrier (Bisnoyatyj^] 
Kogan & Ruzmaikin 1976), and in such a state the inflow must 
somehow accrete through the accumulated magnetosphere. 

The inflow strongly interacts with the magnetospheric barrier 



at the "magnetospheric radius" r m , where the force of gravity by 
the accreting mass equals the magnetic force by the barrier. r,„ is 
readily estimated in a similar way for both neutron stars (NSs) 
( [Lamb et al.|1973|[Elsn"er & Lamb|19T7l|Scharlemann|1978> and 
BHs |Bisnovatyi-Kogan & Ruzmaikin 1974; Nar ayan et al.|2003] >. 
The disk surface density X inside r m is X = M/(2nreVff), where 
M is the mass accretion rate, e = 10e_i is the ratio of mass ad- 
vection velocity to free-fall velocity («ff), where e ~ 0. 1 in our 
simulations discussed la ter. Also, let M = Mh (r/r g )" for horizon 

(1999| . The con- 



accretion rate M H , as in Blandford & Begelman 



dition for magnetic support against disk gravity in the equatorial 
plane is GMX/r 2 ~ 2B r B-J(An) ~ B 2 z /(2n) for B r ~ B z . Then, 
B- ~ 10 5 e: i 1/2 ra 8 1/2 m 1/2 (r/r g )- 5/4+ " /2 G, where r g = GM/c 2 , c is 
the speed of light, G is the gravitational constant, M is the BH 
mass, m g = M/(10 8 M o ), m = Ml Mem, and M Edi = L Edd /(0.1c 2 ) « 
1.4 x 10 25 m 8 g s" 1 is the Eddington mass accretion rate. If e is in- 
dependent of r, integrating B z over r gives an estimate of O within 
the "magnetospheric radius" r m . Inverting <I>(r„,) gives 



112000 



3 n 

4 + 2 



0.1pc 2 Gj 



(1) 

such that r m varies weakly with parameters for higher n and for 
n = (0, 1,2) the coefficient is ~ ( 10 5 , 10 3 , 10 2 )r g , respectively. As 
another measure of how important the magnetic flux is, we also 
consider the limit that all the surrounding flux (O) reaches the BH. 
Then, the horizon's dimensionless magnetic flux would be 



5(r 2 cM H ) 1/2 



10 V 



-3/2 . -1/2 



0.1pc 2 Gj 



(2) 



( |Gamm ie 1999, Penna et al. 2010), which measures the mass- 
loading of the magnetic field lines. An MRI-driven MHD turbulent 
disk has T H < 1 for integrals within the heavy disk inflow ( |Gammie] 
|1999|[Penna et al.|2010| l. If T H <; 1 for integrals over some portion 
of the horizon (e.g. polar regions or heavy disk inflow), then this 
indicates the formation of a force-free magnetosphere and the BZ 
effect can be activated there (Komissarov & Barko vp009^ . 

The quantities r£ , r^ =1 , and Th can be estimated for vari- 
ous systems (using O estimated in a preceding paragraph) to check 
whether a magnetosphere could dominate the flow dynamics. For 
M87 with O ~ 0.1pc 2 G, M « 6. 4 x 10 9 M Q , mass accretion 
rate m H ~ 10~ 4 l |Dexter et al]|201l| , one obtains f'f ~ 10%, 
r" ' " ' ' ' " 

' III 

M 



W 2 r g , and T H ~ W 1 . For SgrA* with $ ~ 0.1pc 2 G 
4.5xl0 6 M o , mass accretion rate m H ~ 1.5x 10 6 ( [Dexter et aT7| 



2010), one obtains r"~ 



10 n r o 



Wr g , and Y H ~ 10 9 . A 
single O-type star feeding SgrA* would give r m ~ r g for n = 0, 1, 2 
and T H ~ 2. For GRS1915+105 with O ~ 10^ I2 pc 2 G, M ~ 15M , 
0.7 ( [Greiner et al.|200l) 



10 3 ^ 



, and T H 
10" 9 pc 2 G, 



Trr 

M 



-14,. 



=1 



For cosmological long-duration GR Bs with 

3M G and M H ~ 0.1M o s" 1 " 
l 



(Woosley et al. 



1993), one obtains r^ =u < r g , r^ 1 ~ r g , and T H ~ 1. One obtains 
only smaller r„, and T H for short-duration GRBs as modelled by 
NS-NS or BH-NS collisions that reach higher M at similar O. For 
GRBs, at late time M H drops, after which r m and T H can be much 
higher (Proga & Zhang 2006). These estimates show that one must 
consider how plasma accretes through a magnetospheric barrier. 
For SgrA* (where n ~ 1 is plausible), even a millionth of O is 
sufficient to lead to magnetospheric accretion near the BH. 

Such flows with accumulated magnetic flux, called "magnet- 
ically choked accretion flows" (MCAFs) by us for reasons based 
upon our simulation results described later, have been considered 
theoretically (Znajek 1976, Bisnovatyi-Kogan & Ruzmaikin 1976; 



© 2012 RAS, MNRAS 000,[T|r35] 



Narayan et al. 2003 1, via pseudo-Newtonian MHD simulations 



( Igumenshchev & Narayan 2002, Igumenshchev et al. 2003; Pen 
|et al.||2003||Igumenshchev||2008[|Pang et al.||2011b, and via gen- 
eral relativistic (GR) MHD (GRMHD) simulations |Tchekhovskoy| 
|et al.|20lT), MCAFs are qualitatively related to flows with strong 



magnetic field (|Shibata et al. 


1990; Mineshige et al.|1995||Machida 


|etal.|2006||Oda etal.|2007, 


Fragile & Meier|2009|. 



MCAFs effectively accrete through a magnetic flux barrier via 
magnetic interchange type modes (Kruskal-Schwarzschild or mag- 



netic Rayleigh-Taylor instabilities) (Wang 1984; Kaisig et al. 


1992 


Lovelace ct al. 


1994 Spruit et al. 1995, Lubow & Spruit 


1995 


Wang 1996; Ste 


hie & Spruit 2001 ; Nakamura et al. 2002; Narayan 


et al. 2003; Li & Narayan 2004, Spruit & Uzdensky 2005; 


Stone 


& Gardiner 2007 , Johansen & Levin 20081. Even with flow shear 



that might remove some specific magnetic interchange modes, still 
other non-axisymmetric magnetic modes operate ( Tagger & Pellatj 
[l999||Stehle & Spruit|2001|>. Interchange modes may also affect jet 
formation and propagation (Nakamura et al. 2002). One might ex- 
pect that magnetic flux accumulation could be prevented by Parker- 
type instabilities (Parker 1966, Johansen & Levin 2008). However, 
the Parker instability cannot operate if the magnetic pressure domi- 
nates the gas pressure due to strong magnetic tension (Shib ata et al.| 
|1990| >. The MCAF-type state is also seen in other systems that de- 
velop a magnetosphere, including young-stellar objects, star for- 
mation regions {Cunningham et al.|2012), an d neut ron stars (|Mes 



tel & Strittmatter 1967 Arons & Lea 1976 1980, Buraard et al. 



1983[|Aly & Kuijpers|1990[|Shu et al.|1994||Romanova et al.|2008| 
Kulkarni & Romanova 2008 , Romanova et al. 201 1). Unlike these 



systems that can have a closed dipolar magnetosphere, BHs harbor 
a split-monopolar magnetosphere (Igumenshchev 2008 I. 

Most 3D pseudo-Newtonian MHD simulations (except in 
Igumenshc hev et al.||2003] [igumenshchev 200 8} and mo st 3D 
GRMHD simulations (except in Tchekhov skoy et al.|2011| |2012| 
|Tchekhovskoy & Mc Kinney 20 f2} of accretion disks have only 
considered relatively weak magnetic fields leading to little accu- 
mulated poloidal flux. Such simulations do not reach the MCAF 
state (for GRMHD simulations, e.g., see |De Villiers et al.||2003| 



McKinney & Gammie|[2004l |McKinney||2005| |De Villiers et al. 



2005; 



Hawley & Krolik||2006| 



McKinney|[2006b| 
20071 



Komissarov & 



McKinney 2007; Komissarov et al. 2007; McKinney & Narayan 
2007a|b| fragile et al.|2007||Beckwith et al.|2008a'fTM cKinney & 
Blandford|2009||Noble et al.|2010||Shiokawa et al.|2012| >, typically 



because their initial conditions used relatively small-sized hydro- 
equilibrium tori within which weak magnetic fields are inserted (so 
only little poloidal flux can accumulate, see section 4.3 in |Igumen^| 
|shchev et al.|2003) . ITchekhovsk oy" et al.| ( |2011| > also used an ini- 
tial torus, but the torus is radially extended and they ensured the 
initial conditions make available a sufficient amount of magnetic 
flux to reach flux saturation near the BH. Note that 3D is required 
for simulations to avoid decaying turbulence (Cowling 1934b and 
to properly resolve a MCAF. For example, the suspended inflow 
seen in 2D axisymmetric MHD simulations (McKinney & Gam- 



|mie|2004| [Proga & Zhang|[2006) |Komissarov & Barkov|2009) is 
due to axisymmetry Pgumenshchev|2 008 2009). Other 2D MHD 
simulations have constant poloidal flux due to using a purely radial 
field ( Proga & Begelman 2003 ), which gives results similar to torus 
simulations with limited poloidal flux. 

Also, MHD simulations have typically studied relatively thin 
disks compared to expected for radiatively inefficient accretion 
flows (RIAFs), as applicable to (e.g.) SgrA*, M87, and some x- 
ray binary states. Quasi-analytical RIAF models include advection- 
dominated accretion flows (ADAFs) ( Ichimaru 1977 , Narayan & Yi 
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l994l|1995b||Abramowicz et al.|1996|[T995||Popham & Gammie 



1998 



et al. 



convection-dominated accretion flows (CDAFs) (Narayan 



2000, Quataert & Gruzinov 2000), and advection-dominated 



inflow-outflow solutions (ADIOSs) (Blandford & Begelma nfl999| 
Begelman 2011). These models suggest that RIAFs should have 
disk height (H) to cylindrical radius (R) ratio of \H/R\ ~ 0.5-0.9. 
Simulations have studied \H/R\ ~ 0.05-0.1 ( |Shafee et aLl|2008l 



Reynolds & Mil ler|2009| |Noble et al.|2009||201 0; Sorath iaeFaT 
2010[ |Beckwith et al ||2011|>, \H/R\ ~ 0.1-0.15 dHawley & Kro- 
|lik|2001||De Villiers et al.|2003||Beckwith et al.|2008b|a| |2009| >7 
\H/R\ ~ 0.2 (|Hawley & Balbus|2002)|Machida et al.|2000[|Machida 



& Matsumoto||2003| |M cKinney & Gammie 2004; Fragil e et al. 
2007[|McKinney & Blandford||2009) , \H/R\ ~ 0.3-0.4 PgumerT 
shchev et al.|2003[|Penna et aLgOlok \H/R\ ~ 0.6 (Run F type in 



Stone & Pringle 2001 ; Hawley et al. 2001 1, and rarely \H/R\ 



1 

i Igumenshchev & Narayan|2002||Pen et al.|2003||Pang et al.|201 1| ). 

In this work, we use fully 3D GRMHD simulations to study 
radially extended and thick (including \H/R\ ~ 1) radiatively ineffi- 
cient BH accretion flows with poloidal and toroidal magnetic field 
geometries and various a/M. Our initially poloidally-dominated 
field models are designed so that near the horizon the poloidal 
magnetic flux reaches a saturation point independent of the initial 
poloidal magnetic flux. Most prior MHD simulations would have 
only required a few times more poloidal magnetic flux to reach this 
natural saturation point. The only other natural limit of poloidal 
magnetic flux is it's negligible. So, for comparison, we also con- 
sider models with initially toroidally-dominated magnetic field. 

The equations solved are presented in ^2] diagnostics are de- 
scribed in ij3] and models are described in ^4] Results for our fidu- 
cial model of a thick disk around a rapidly rotating BH are in ^5] 
Results for various field geometries, BH spins, and disk thicknesses 
are in ^6] We discuss our results in ^7] and conclude in ^8] Some 
numerical method details are given in Appendix [A] 



2 GOVERNING EQUATIONS 

We solve the GRMHD equations for a radiatively inefficient mag- 
netized accretion flow around a rotating black hole defined by the 
Kerr metric in Kerr-Schild coordinates with internal coordinates 



^ = (?,x' 



,(1) r (2) r (3) 



) mapped to the spherical polar coordinates 



rf = (t,r,6,cf>). We write orthonormal vectors as u h contravariant 
(covariant) vectors as m' (m.), and higher-ranked coordinate basis 
tensors with no underbar. We work with Heaviside-Lorentz units, 
often set c = GM = 1, and let the horizon radius be r H . 
Mass conservation gives 



V„(p M") = 0, 



(3) 



where po is the rest-mass density, if is the contravariant 4-velocity, 
and p = p u' is the lab-frame mass density. 
Energy-momentum conservation gives 



v„r? = 0, 



(4) 



where the stress energy tensor 7f includes both matter (MA) and 
electromagnetic (EM) terms: 

r MA '; = (po + u g + Pg )u"u y + p g 5» v , 



rEMf _ ,2u 



b Z U»U y + p„^ - ^b v 



(5) 



The MA term can be decomposed into a particle (PA) term: T PA y = 
Po l l v \f an d an enthalpy (EN) term. The MA term can be reduced to 
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a free thermo-kinetic energy (MAKE) term, which is composed of 
free particle (PAKE) and enthalpy (EN) terms: 



r MA C- Po ^ la, 



(6) 



rPAKE'' 



(u y - rj_Ja)p if, 

T EN " V = (u g+Pg )i/u v + p s S>;„ 

such that T MAKE 1 = T FAKE " V + T w ^,. Here, u g is the internal energy 
density and p g = (T — l)« g is the ideal gas pressure with adiabatic 
index T = 4/3 (r = 5/3 may lead to somewhat different results 
; |McKinney & Gammie(2004) |Mignone & McKinney|2007| >. The 
contravariant fluid-frame magnetic 4-field is given by If, which is 
related to the lab-frame 3-field via If = ff'l^/u' where ]f v = ifu y + 
6y is a projection tensor, and dy is the Kronecker delta function. 
The magnetic energy density (ut) and pressure (pt) are uj, = pi, = 
lfb/2 = b 2 12. The total pressure is p toi = p g + p b , and plasma 
P = Pg/pb- The 4- velocity of a zero angular momentum observer 
(ZAMO) is i] = {-a, 0, 0, 0) where a = 1/ -y^g" is the lapse. The 
4-velocity relative to this ZAMO is if = if-yrf where y = -u a T] . 
For any 3-vector (e.g. 6'), the "quasi-orthonormal" vector is S, = 
ff y[g[j computed in spherical polar coordinates. 

Magnetic flux conservation is given by the induction equation 

d, ( V=gfi') = -9j[ V - B V)] , (7) 

where g = Det(g /JV ) is the metric's determinant, and the lab-frame 
3-velocity is «' = m'/h'' No explicit viscosity or resistivity are in- 
cluded, but we use the energy conserving HARM scheme so all 
dissipation is captured (Gammie et al. 2003, McKinney 2006a). 

The energy-momentum conservation equations are only mod- 
ified due to so-called numerical density floors that keep the numer- 
ical code stable as described in detail in Appendix|A] The injected 
densities are tracked and removed from all calculations. 



3 DIAGNOSTICS 

Diagnostics are computed from snapshots produced every ~ 2r g /c. 
For quantities Q, averages over space {{Q)) and time ([<2D are per- 
formed directly on Q (e.g. on rather than on any intermediate 
values). Any flux ratio vs. time with numerator F N and denomi- 
nator F D (F D often being mass or magnetic flux) is computed as 
R(t) = (F N (f))/[{F D y] t . Time-averages are then computed as [/?],. 

3.1 Fluxes and Averages vs. Radius 

For flux density F d , the flux integral is 

F(r) = J dA 23 F d 
where CIA23 = ^f^gd^d^ (dA g $ is the spherical polar version) 



(8) 



For example, F d = po^ (1> gives F = M, the rest-mass accretion rate. 
For weight w, the average of Q is 

Q w (r) = <e>„. = Vtt^ — . (9) 



All 6, (f> angles are integrated over. 



3.2 Fluxes and Averages vs. 9 

The flux angular distribution, at any given radius, is 

Xx/2-\6-it/2\ r-n 
dA^F d + dA^F d , (10) 

=0 Je'=i/2+|fl-ff/2| 



which just integrates up from both poles towards the equator, is 
symmetric about the equator, and gives the total flux value at 6 = 
7r/2. The average of Q vs. 6 using weight w is given by 



QJ6) 



J dAg^wQ 



(11) 



All (^-angles are integrated over. 



3.3 Disk Thickness Measurements 

The disk's geometric half-angular thickness is given by 



«<»-*af. 



(12) 



where we integrate over all for each r,(f>, and 6 = n/2 + 
{(8 — n/2)) p is also integrated over all 6 for each r, tp, and the fi- 
nal ff'(r) is from ^-averaging with no additional weight or -\/-g 
factor. This way of forming ff'ir) works for slightly tilted thin disks 
or disordered thick disks. For a Gaussian distribution in density, 
this satisfies p/(p[9 = 0]) ~ exp(-0 2 /(2(0 rf ) 2 )). For sound speed 
c s = *\/Tp g /(p + u g + p g ), the thermal half-angular thickness is 



= arctan 



(c s ) w 



(13) 



where v 2 ot = 
plerian (i.e. 



+ Vg. For a thin hydrostatic non-relativistic Ke- 
= \v K \ with v K a R/(a + R 3 ' 2 )) Gaussian disk, 
(f = c s /w ro t for c s and v mt at t he disk plane. Also, ff p ~ 0.93 0*' for 
T = 4/3. Ram pressure forces (jBeskin & Tc hekhovskoy 2005 1 and 
magnetic forces l Cao|20lT 1 can cause 1 <K 0'. See Tables [T] [3] 
Note that ADAFs have ff > 1 jNarayan & Yi|1994f . 



3.4 BH, Disk, Jet, Magnetized Wind, and Entire Wind 

Many quantities (Q) vs. r or vs. 6 or vs. are considered for various 
weights and conditions. We define the superscript "f" (full flow) 
case as applies for weight w = 1 with no conditions, "fdc" (full 
flow except avoids highly magnetized jet where numerical floors 
are activated), "dc" (disk plus corona but no jet) case as applies 
for w = 1 with condition b 2 /p < 1, "dcden" (density- weighted 
average) with w = p and no conditions, "ff 1 " (within 1 disk half- 
angular thickness) case with w = 1 and condition of \0 - 6q\ < ff 1 , 
"eq" (within 3 cells around the equator) case with w = 1, and "jet" 
or "j" case (jet only) with w = 1 and the condition that density 
floors are activated (see Appendix[A|. For quantities vs. 6 or vs. <f>, 
we radially average within ±0.1 r at radius r. 

Fluxes, described in the next section, have integrals computed 
for a variety of (somewhat arbitrary) conditions. The subscript 
"BH" or "H" is for all angles on the horizon. The subscript "j" or 
"jet" is for the "jet" with condition b 2 /p > 1. When the jet is mea- 
sured at a single radius, we use r = 50r g (except the MB09Q model 
that uses r = 30r g due to its limited radial range). The subscript 
"mw" is for the "magnetized wind" with conditions b 2 /p < 1 and 
p < 2 for all fluxes, except for the rest-mass flux that also has 
-(po + u g + P g )u,lpo > 1 (i- e - thermo-kinetically unbound). The 
"w" or "wind" subscript is for the "entire wind" with the condi- 
tion of b 2 /po < 1 that includes all of the flow except the jet. The 
subscript is "in" ("out") for the condition u r < (u r > 0). 

We also compute (as shown in Table |3j the number of grid 
cells that cover the disk half-angular thickness (ff 1 ) computed by 
Eq. |12| at the horizon (ff^) and other radii r (6f /r ). For compari- 
son, we also compute the thermal half-angular thickness (ff = ff ). 
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Also computed are the disk-corona interface angular location at 
a given radius r (denoted 6rj r ) defined by where B = \ and the 
corona-jet interface angular location at a given radii r (denoted 
(f' lr ) defined by where b 2 /p = 1. In practice, these interface lo- 
cations are defined similarly to Eq. (12) , except 9q = n/2 and 
weight w = «'(po + u g + p g + b 2 ) are chosen with the following 
conditions. The disk-corona interface calculation uses the condition 
1/2 < 8 < 1, unless that condition is not met by any grid cells at 
that radius - in which case the condition 1/10 < B < 1 is used. The 
corona-jet interface calculation uses the condition 1 < b 2 /p < 2, 
unless (very rarely) that condition is not met by any grid cells at 
that radius - in which case the condition 30 > b 2 /po > 1 is used. 



3.5 Fluxes of Mass, Energy, and Angular Momentum 

The rest-mass flux, specific energy flux, and specific angular mo- 
mentum flux are respectively given by 



J 



M 

E 

J 



J p u r dA e , 



and are computed in Tables[3J|6] 

The net flow efficiency is given by 

E-M E EM (r) + £* 



im, 



[M HJ( 



(14) 
(15) 
(16) 

(17) 



Positive values correspond to an extraction of positive energy from 
the system at some radius. These ^'s are computed in Tables[5][6] 
The BH's dimensionless spin-up parameter is 



d(a/M) M 



-1-2 — (l-f7), 



(18) 



(computed in Table |7|. All 9 and <f> angles are integrated over. The 
BH is in "spin equilibrium" for s = (Ga mmie et al.|2 004l. 



3.6 Magnetic Flux 

The radial magnetic flux vs. 9 at any radius is 



V r (r,8)= dA e4 ,B r 



(19) 



The signed value of the maximum absolute value over all 6 angles 
(smaxa e ) of the magnetic flux is 



*P t (r) = smaxao 1 ?,-, 



(20) 



and = l^r = ^h) > s the horizon's magnetic flux. The half- 
hemisphere horizon flux is 



T H = *F r (r = r H ), 



(21) 



as integrated from 6 = n/2 to n (negative compared to the integral 
from 9 = to n/2). The 9 magnetic flux vs. radius at angle 9 is 



y e {r,9)= f V = ^ (1) ^ <3> S- 2> . 



(22) 



where the vertical magnetic flux threading the equator is 

y eq (r) = y e (r,e = n/2). (23) 



The total magnetic flux along the equator is 
f(r) = T H + 4% (r). 



(24) 



For all forms of T*, all 0-angles are integrated over. 

The magnetic flux can be normalized in various ways (as com- 
puted in Table |9j. Normalization by the initial flux at r gives 
v F(r)/ v P(r ). One type of field geometry we will use has multiple 
field loops of alternating polarity as a function of radius. So an- 
other normalization is by the initial i-th extrema vs. radius, which 
gives f/f, that picks up the extrema in the magnetic flux over each 
field loop. Normalization by the initial value of an extrema gives 
y/^iit = 0). We also need to form a measure that indicates how 
much flux is available to the BH. So we consider the normalization 
by the flux in the disk that is immediately available to the horizon 
of the same polarity. This measure is given by Yh/Y,,, where *¥„ is 
the value where T(r) goes through its first extremum of the same 
sign of magnetic flux (i.e. out to the radius with the same polarity of 
dipolar-like field) as on the horizon. If the horizon value is itself an 
extremum, then Th/T,, = 1 implying that the region immediately 
beyond the horizon only has opposite polarity field. 

The absolute magnetic flux (O) is computed similarly to X V, 
except one 1) inserts absolute values around the field (e.g. B' and 
S" in the integrals); 2) puts absolute values around the integral ; 
and 3) divides by 2 so that a dipolar field has ITJ = O. For exam- 
ple, <S r (r,9) = (1/2) |J dA^\ff\\. The quantity O/'F, (computed in 
Table [5] and which is the only flux ratio directly time-averaged as 
[O/T*,],) is roughly the vector spherical harmonic multipole / of the 
(^-component of the magnetic vector potential: 



-^4 



t 

Jfl' = 



as integrated over all </>. For example, for / 
10/^,1 = 1, 2, 2.6, 3.5, 4, 5.6, 5.7, and 6.7. 

The Gammie ( 1999) model normalization gives 

Y * 0.7 ' . 



(25) 

1 . . . 8) one gets 



(26) 



which accounts for O r being in Heaviside-Lorentz units ( |Penna| 
|et al.| |2010). Compared to Gaussian units version of (/>h = 



4W- 



Tchekhovskoy et al. 



(2011 



a O.20H- 
and T w 



T H and T ; are normalized by M H , T in by M in , and T ra 
respectively by M mw and M w . T is computed in Table]?] 

The field line rotation frequency with respect to the BH spin 
(z) axis is computed various ways. We consider fip = F tr /F rl/ „ fip = 

Fa/F^, fi F = \v*\ + S ign[u r ](v p /B„)\Bt\ with v p = ^jv 2 + v] and 
B p = tJb 2 + B], and 



VyBy + VyBf, 

m+Bi 



(27) 



We also consider £2 F = [|F, e |],/[|F e(4 |],. These flp are normalized by 
the BH rotation angular frequency Q H = a/(2Mr li ). 



3.7 Inflow Equilibrium and a Viscosity 

Inflow equilibrium is defined as when the flow is in a complete 
quasi-steady-state and the accretion fluxes are constant (apart from 
noise) vs. radius and time. The inflow equilibrium timescale is 



: N 



(28) 
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for N inflow times from r = r K and r, = \2r„ to focus on the 
more self-similar flow. t, is used in Table 



1, andr? 



10 



uses r,, with N = 3. 



Viscous theory gives a GR a-viscosity estimate for v r of 

v vhc Ga(6 d ) 2 \v m \ ( |Page & Thome||l974[|Penna et al.||2010| >, 

with GR correction G (< 1.5 for r > 58r g ) and (not the lapse) 



(29) 



<*PA 



«M2 



ffeff 



p 6u r {8u.Jg^) 



Plot 



Ptol 



P6 



p g )&u r {6u 4 



"vise /a 

a e ff(|Vro 



.iI/IwkI), and (small) a EN a (i^ + 
and a m ~ tP-du^du^ yfg^)/p tat - Here, 
Su is the deviation of the velocity from its average (taken over 
all <f> and over the time-averaging period). The a (e.g. in Table [5| 
is averaged as follows. The numerator and denominator are 
separately volume averaged in 9, <f> for each r. Weight w = 1 with 
condition b 2 /po < 1 gives a a for the disk+corona, while w = p 
gives a b for the heavy disk. Notice a M2 = <*mag/(l +Anag) for some 
P denoted j8 mag , and sm(26 b ) = a raag for tilt angle 9 b (jSorathia et al 



2011). These o-'s are accurate for \v\ « cas true for r > 2r g in 
our models, while a cS is accurate far outside the inner-most stable 
circular orbit (ISCO). 

3.8 Modes and Correlation Lengths 

The flow structure is studied via the discrete Fourier transform of 
dq (related to quantity Q) along x = r, 9, if> giving amplitude a p for 
p = n,l, m, respectively. The averaged amplitude is 



p \) = (\r p (dq)\)= f ^dq 

J not x . n 



(30) 



computed at r = r H , 4r g , Sr g , 30r g . The x is one of r, 9, cf> and "not 
x" are others (e.g. 9, <f> for x = r). The dq is (generally) a function 
of x on a uniform grid indexed by A: of A' cells that span: Sr equal 
to 0.75r around r for x = r, n for x = 9, and 2n for x = <f>. The N is 
chosen so all structure from the original grid is resolved, while the 
span covered allows many modes to be resolved. 

For all x, dq = -s/^gdx^dx^dx^SQ/qN. For x = r,9, we let 
9» = Lot x ^Fgdx m dx {2) dx & ([Ql), ([QD as the time-0 averaged 
Q, and SQ = Q - ([Q] t )- Using dq removes gradients with r, 9 
so the Fourier transform acts on something closer to periodic with 
constant amplitude (see also Beckwith et al. 2011). For x = cf>, 
we let q N = 1 and SQ = Q because the equations of motion are 
0-ignorable. For x = 9,<f>, the radial integral is computed within 
±0.1 r. For x = r, 9, the <p integral is over all 2n. For x = r,cf>, 
the 9 integral is over all n. For all x cases, the 9 range of values 
uses the "fdc" or "jet" conditions (respectively called "Disk" and 
"Jet" in sections ["5. 7|6,7|6,l0| l, where these conditional regions are 
defined via ^-averaged quantities at each time. Notice we average 
the mode's absolute amplitude, because the amplitude of (SQ) de- 
resolves power (e.g. m = 1 out of phase at different 9 gives (SQ) — > 
and a„, — > 0) and is found to underestimate small-scale structure. 

We also compute the correlation length: A y col = x cor - Xq, 
where x = for x = 9, if) and x is the inner radius of the above 



given radial span for x = r, where rc cor = Sr/A riCOr , l cm = 7r//l fliCor , 
and m cal = (2n)/A,p cm . The Wiener-Khinchin theorem for the auto- 
correlation gives 



exp(-l) : 



T-J X0 [\(a p>0 \) 2 ] ' 



(31) 



where T 1 [(|a ;)> ol) 2 ] is the inverse discrete Fourier transform of 
(|a ; ,|) 2 but with (<7o) reset to (i.e. mean value is excluded). 

3.9 Suppression of the MRI 

The MRI is a linear instability with fastest growing wavelength of 

\Vx,k\ 



« 2n 



|OJ' 



(32) 



for x = 9, (f>, where \v x _ A | = Jb^WJe is the x-directed Alfven speed, 
e = b 2 +po+Ug+p g , and rQ. lot = v mt . Ausa is accurate for £l m , oc r~ 5/2 
to . Q rot , Va are separately angle-volume-averaged at each r, t. 

The MRI suppression factor corresponds to the number of 
MRI wavelengths across the full disk: 

2r9 d 

SdMsi - ~j • (33) 

Wavelengths A < 0.5/1 9 mri are stable, so the linear MRI is sup- 
pressed for Srf.MRi < 1/2 when no unstable wavelengths fit within 
the full disk jBalbus & Hawley|[i998| jPessah & PsaTtisl[2005) ). 
Sdum Srf.weatMRi) uses averaging weight w = (b 2 p) 1 ' 2 (or 
w = p), condition j8 > 1, and excludes regions where density floors 
are activated. Weight w = (b 2 p) [/2 is preferred, because much mass 
flows in current sheets where the magnetic field vanishes and yet 
the MRI is irrelevant. When computing the averaged S^mri, f a and 
|fl rot | are separately 9, <p- volume- averaged within ±0.2r for each t, r. 
The averaged S^mri is at most 30% smaller than 5,;, W e a k,MRi- 

S rf ,f=o,MRi,{i,oi (in Table[T} gives 5 rf , M Ri at t = at r = n = 30r g 
and r = r„ = 50r g averaged within ±0.2r, except models MB09D/Q 
use r, = 10r g and r = 15rpdue to their disk's limited radial extent. 
Time-averages (see Table [8| are obtained for r = rf 



(see Table 10 1 averaged within ±0.2r. Table[8]also gives , 
r Si/MRi=i/2' within which the linear MRI is suppressed. 



4 PHYSICAL AND NUMERICAL MODELS 

This section describes our models with parameters shown in Ta- 
bles [T| [2] The model names are in the form AxByNz, where x is 
the approximate value of the BH spin, y identifies the field geom- 
etry (p=poloidal, f=flipping poloidal, t=toroidal), and z identifies 
the normalization of the magnetic field. For instance, our fidu- 
cial model A0.94BfN40 has a spinning BH (a/M = 0.9375), a 
poloidal field that flips polarity with radius, and /? m j n « 40 (i.e. 
smallest value of /? is 40). The label c? (with number ?) is ap- 
pended for convergence tests, r is appended to the model name 
if it is another realization of an identical model, and HR is ap- 
pended if it is a high-resolution continuation of some model. 2D 
axisymmetric models are marked with *. Models from McKinney 
|& Blandford] ( [2009) are denoted MB09D and MB09Q for their 
dipolar and large-scale quadrupolar models, respectively. The re- 
maining TMN1 1 models are similar to Tchekhovskoy et al.| l |20PT| >. 
Primary models (fiducial thick poloidal: A0.94BfN40, thick ret- 
rograde poloidal: A-0.94BfN40HR, thick toroidal: 0.94BtN10HR, 
thinner poloidal: A0.99N100) have bold font labels. 

Table [2] gives the time-averaging period (from Tf to T" f in 
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r g /c). Models have T t = 0, except A0.94BtN10HR starts as a 
super-sampled A0.94BtN10 at T, = 57400r ? /c, A-0.94BfN40HR 
starts as A0.94BfN40 at T t = 7995r g /c with a/M = -0.9375, 
and A0.99N100 doubles 2V # (to shown A^) at t « 14675^/c. For 
A0.94BpN100, T%V f = 14000,27048 and Tf, T" = 8000, 13000 
give similar results (validating Tj = 13000 for similar models). 

4.1 Physical Models 

This study considers BH accretion disk systems. Our "thick 
disk" models have initial geometric half-angular thickness 0* ~ 
0.6, a range of BH spins (a/M = - 0.9375, 
0.5, 0, 0.5, 0.9375), a range of field geometries (constant polar- 
ity poloidally-dominated, flipping polarity poloidally-dominated, 
and toroidally-dominated), and a range of initial disk magnetic 
field strengths (J3 mm = 10, 30, 40, 100). Our "thinner disk" 
TNM11 models have initial ff' ~ 0.2, a/M = - 0.9, - 0.5, - 
0.2, 0, 0.1, 0.2, 0.5, 0.9, 0.99, a poloidally-dominated field ge- 
ometry, and /3 mi „ = 25, 50, 100, 200 for a/M = 0.9. MB09 models 
are like most prior MHD simulations with limited poloidal flux. 

The initial mass for the thick disk models is an isentropic 
hydro-equilibrium torus ( jFishbone & Moncrief] [T976] |Gammie| 
|et al.|2 003 ) with inner edge at r m = 10r g and pressure maximum 
at r max = 100r s . The torus is marginally unbound by tens of per- 
cent, as similar to ADAFs that we want to model ( |Narayan & Yi| 
|1995a[ l. The magnetic field inserted (described later) makes negli- 
gible changes to the torus' boundedness. Table[3]shows the disk's 
geometric half-angular thickness (ff 1 ) at r max . We set p = p max = 1 
at the maximum rest-mass density. To seed the MRI, u g is perturbed 
by a factor 1 + F R (E — 0.5), where F R = 0. 1 and £ is a random num- 
ber from to 1. The torus is surrounded by an atmosphere with 
p = 10r*(r/r g )- 2 , u g = \Q- b (r/r s y sl2 , ff = 0, and B> = 0. 

We consider an initial poloidal field geometry to seek poloidal 
magnetic flux saturation near the BH. Field polarity flips are in- 
serted by modulating the poloidal polarity vs. radius in order to 
generate multiple loops of alternating polarity for studying mag- 
netic field inversion/annihilation. For this poloidal field geometry, 
the ^-component of the magnetic vector potential is 

A* <* /1/2. (34) 
/, = M'Vsin(0)r, 

*? ^£,max Jo 

f 2 = sin (log (r/S)/:T), 

where f\ has p = 1 and v = 2, q has f c = 0.2, u gtmslx is the maximum 
u g , q = is set if q < 0, and f 2 has 5 = 0.5r m and T = 0.28 for the 
flipping field and f 2 = 1 for the non-flipping field. 

We also consider an initial toroidal field geometry, as the 
limit of negligible coherent poloidal magnetic flux, where the 0- 
component of the magnetic vector potential is 

A g oc r 2 (u g /u m . ix - f c ), (35) 

with f c = 0.2. If Ag < 0, then A fi = is set. Then V^S^ = 
Ag r is computed. The random perturbations of u g also lead to a 
small radial field via = - A^, which corresponds to radial 

wiggles in the toroidal field. Within r ~ 100r s there is about 10 
times less energy in this radial field compared to the toroidal field. 
Very small truncation-level B fl is also present. 

Table[T]shows each model's initial field, marked as "Poloidal" 
for the non-flipping poloidal field, "PoloidalFlip" for the flipping 
poloidal field, and "Toroidal" for the toroidal field. The mod- 
els from |McKinney & Blandford] p009"l > use the "PoloidalOld" 



poloidal field geometry (i.e. oc (p - p cut ) and, e.g., p cut ~ 
0.25p m;lx = 0.25) or use the "LSQuad" large-scale quadrupolar field 
geometry. The models marked as "Poloidal2" use the magnetic field 
geometry described in Tchekhovskoy et al. (201 1 1. 

The magnetic field strength is set via the plasma f} = 
PglPb ~ (2/T)(c s /v cl f where v 2 = b 2 /(p + u g + p g + b 2 ) gives 
the Alfven speed v a . Our thick disk models have yS m j n , the smallest 
value of fi (within the resolved disk region, e.g., r ~ 1000r s <c R out ) 
of /3min ~ 10 to 200. An alternative measure is y6 rat - f-maxes = 
Ps,max/Pfc,max> where p s _ max is the maximum thermal pressure on the 
domain and pi. max is the maximum magnetic pressure on the do- 
main. Another alternative is /3 m t-of-av g = P g .avg/Pi,.av g = (p g )/(Pb)- 
These /? (see Table [TJ are computed with condition b 2 /p < 1. For 
poloidal field models, our choices for fi ensure that S^mri > 1 so 
the MRI operates, while we push close to 5,;_mri ~ 1 as reached 
even in toroidal field models. 

4.2 Numerical Models 

The uniform spatial coordinates x (,) have resolution A',, x N e X 
active grid cells and 4 boundary cells for each of the 6 boundaries in 
3D. The radial grid of N r cells spans from R in to i? out with mapping 

r(x m ) = R + exp f[x m ] (36) 

where Rq = is chosen in this paper. For x m < x\, IC ak 

/[x (1) ] = n x m , (37) 

where Xb rea k = log(rb rea k _ Ro)/ n o (with no = 1), and otherwise 

f[x w ] = n x (l > + c 2 (x m - x brcak r, (38) 

where c 2 = 1, and n 2 = 10. The x (1) grid ranges from x^ 1 ' = 
(log(i? ta -2?o))/no tox^ , which is x ( f [) = (log(R out -R ))/n ifR out < 
''break and otherwise determined iteratively from R 0M = r[xr?]. The 
value of R[ n is chosen so that there are 6 active grid cells inside 
the outer horizon, while R iu is outside the inner horizon. So the 
boundary cells only connect to stencils (each ±4 cells) that are in- 
side the horizon, which avoids causal connection between the inner 
boundary and the flow outside the horizon. For models where no 
persistent jet is launched, we set R out = 10 3 r e and r brea k » R mt . 
For models where a jet is launched, we set R out = 26000r s and 
'"break = 5 x 10V S . The radius r break is where the grid changes from 
exponential to hyperexponential, which allows the grid to focus on 
the dynamics at small radii while avoiding numerical reflections off 
the outer grid. Radial boundaries use absorbing conditions. 
The 0-grid of N e cells spans from to n with mapping 

6(x a) ) = 2 + W(e, - 2 ) (39) 

where x <2) ranges from to 1 (i.e. no polar cut-out ; but see Ap- 
pendix[A|. The first grid mapping function is given by 



0, = 




T 2 ~- 


= (n/2)(l + arctan (h 2 (x {m2) - (1/2)))/ arctan (h 2 /2)), 


h 2 ~- 


-- h+((r-r sji )lr 0ji fi\ 


To -- 


-- nx\ 2) + ((1 - h ) /2) sin (2^x <2) ), 


ho 


_ 2— Q (r/ri .) _ "j5"'' 2 l | ' l ' I '" rdIll ''''«r'>i''i/ll 


s 2 -- 


= (1/2) -(1/7T) arctan ((r-r s )/r ), 


So ~- 


-- (1/2) + (1/ W ) arctan ((r-r s )/r ), 



where r s = 40 and ro = 20. For h 2 , we set A3 = 0.3, rofl = 20, 
r s p = 0, and nji = 1 so the jet is resolved with grid lines following 
9j oc r~"P . For h , we set r,^ = 2.8, np = 1, r j = 15, r s t = 40, and 
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Table 1. Physical Model Parameters 



ModelName 


a/M 


FieldType 




/^rat-of— avg 


^rat-of-max 


< m 


6>;. mai 


5rf,(=0,MRI,(i,ol 


Tf 


A A II II? CXI Uk 


U.93 ij 


PoloidalFlip 


40 


230 


40 


0.61 


i.i 


1,0, 2.0 


26548 


AU.y4iilJN 1 UUC 1 


0.9375 


PoloidalFlip 


120 


ion 

4yu 


1 30 


O.o 1 


i.i 


"t A A A 

2.4, 4.4 


1 3000 


AU.y4±3tiN 100C2 


U.V3 /j 


PoloidalFlip 


1 30 


4MJ 




0.0 1 


■ 


1 1 A A 
2.3, 4.4 


1 3U00 


A0.y4BIJN 100C3 


0.9375 


PoloidalFlip 


i in 
IzU 


4VU 


1 30 


O.o 1 




2.3, 4.5 


1 3000 


A0.94BfN100c4 


0.9375 


PoloidalFlip 


130 


490 


130 


0.61 


u 


2.3, 4.4 


13000 


A0.94BfN40c5* 


0.9375 


PoloidalFlip 


40 


230 


43 


0.61 


i.i 


1.6, 2.6 


13000 


A-0.94BfN40HR 


-0.9375 


PoloidalFlip 


40 


230 


40 


0.61 


i.i 


1.6, 2.6 


18416 


a r\ O/iDfMin 
A-U,y4oiJNjU 


-0.9375 


PoloidalFlip 


3U 


1 20 


3 1 


U.o 


i.i 


1.2, 2.2 


1 3000 


A-U.jolJNiU 


-0.5 


PoloidalFlip 


34 


i it) 
13U 


35 


0.61 


i.i 


1.3, 2.2 


1 3000 


AU.UdiJN 10 





PoloidalFlip 


1 


42 


1 1 


0.61 


i.i 


0.68, 1.2 


13832 


A0.5B1JN30 


0.5 


PoloidalFlip 


33 


120 


36 


0.61 


i.i 


1.3, 2.2 


1 3000 


A0.y4BilN30 


0.9375 


Foloidalrlip 


30 


120 


32 


0.61 




1.2, 2.2 


13000 


A0.94BfN30r 


0.9375 


PoloidalFlip 


32 


120 


33 


0.61 


u 


1.2, 2.2 


13000 


AU.y4BpN 100 


0.9375 


Foloidal 


93 


240 


1 10 


0.61 




1.8, 1.8 


27048 


A-0.94BtN10 


-0.9375 


Toroidal 


10 


1600 


25 


0.6 


1.1 


370, 350 


93944 


A-0.5BtJN 10 


-0.5 


Toroidal 


10 


980 


17 


0.61 


1.1 


300, 260 


129392 


AU.uBtJN 10 





loroidal 


10 


1200 


28 


0.61 


1.1 


Tin T f\r\ 

320, 300 


96796 


A0.5BtN 10 


0.5 


Toroidal 


to 


1000 


17 


0.61 


1.1 


290, 250 


135572 


AU.y4BtJN 1U 


0.9375 


Toroidal 


10 


1300 


20 


U.ol 


1.1 


3UU, ly\J 


93648 


A A I I 1 A f 1 1 > 

AU.94r5tN luHK 


0.9375 


Toroidal 


10 


1300 


20 


0.61 


1.1 


inn or\r\ 

300, 290 


80184 


MBOyD 


0.92 


FoloidalOld 


14 


530 


100 


0.18 


0.24 


3.4, 16 


5662 


JVlBUyv,) 


0.9375 


Loljuaa 


0.0067 


540 


140 


0. 19 


0.25 


41, 60 


jd88 


a r\ r>\! i Art 

A-0.yN100 


-0.9 


Foloidali 


1 00 


150 


220 


0.21 


0.27 


3.6, 4.4 


20060 


A A CAT 1 r\r\ 

A-0.5N100 


-0.5 


Foloidali 


1 00 


140 


210 


0.22 


0.28 


3.5, 4.3 


16335 


A-0.2N100 


-0.2 


Poloidal2 


100 


140 


200 


0.22 


0.29 


3.5, 4.3 


15175 


A0.0N100 





Poloidal2 


100 


140 


210 


0.22 


0.29 


3.6, 4.3 


18585 


A0.1N100 


0.1 


Poloidal2 


100 


150 


210 


0.22 


0.29 


3.6, 4.4 


20000 


A0.2N100 


0.2 


Poloidal2 


100 


140 


210 


0.23 


0.29 


3.6,4.3 


19350 


A0.5N100 


0.5 


Poloidal2 


100 


150 


210 


0.23 


0.29 


3.6, 4.4 


19540 


A0.9N25 


0.9 


Poloidal2 


25 


36 


51 


0.23 


0.3 


1.8, 2.2 


17350 


A0.9N50 


0.9 


Poloidal2 


50 


72 


100 


0.23 


0.3 


2.5, 3.1 


14385 


A0.9N100 


0.9 


Poloidal2 


100 


140 


210 


0.23 


0.3 


3.6, 4.4 


19895 


A0.9N200 


0.9 


Poloidal2 


200 


290 


420 


0.23 


0.3 


5.1, 6.2 


28620 


A0.99N100 


0.99 


Poloidal2 


100 


160 


230 


0.24 


0.3 


3.7, 4.6 


31400 



Qj = 1.3. Wesetx (m2) = x (2) unless x (2) > 1 thenx ( "' 2) = 2-x (2) and 
unless x (2) < then x (m2) = -x <2) . ®i focuses on the disk at small 
radii and the jet at large radii. The second mapping function is 

0, = (n/2)(h e (2x i2) - 1) + (1 - /i fi )(2x (2) - I)"* + 1), (41) 

where n g = 5 and h e = 0.15. ®2 focuses on the thin inflow near the 
horizon in poloidal field models, while it also avoids small <p polar 
cells that would limit the time step. The interpolation factor is 

W = (1/2) + (l/;r)(arctan((r - r s j 2 )/r 0j2 )), (42) 

where r s p. = 5 and roj2 = 2. The polar axis boundary condition is 
transmissive as described in Appendix [A] 

The 0-grid of N$ cells spans from to 2n with mapping 
0(x <3) ) = 27tx (3) . Many of our simulations have x (3) vary from 
to 1 such that A(f> = 2n. This is a fully 3D (no assumed symmetries) 
domain. Periodic boundary conditions are used in the (^-direction. 
Our TNM11 type models use various A0, and the spatial integrals 
are renormalized to account for the full 27r range in <p. 

Table [2] marks the grid as "Exp" if exponential, "HypExp" 



if hyperexponential, "ExpOld" for MB09's exponential grid, and 
"TNM11" for TNM1 1 's hyperexponential grid. 

We choose a resolution N r x Ng x N$ that has a grid aspect 
ratio of 1 : 1 : 1 for most of the inner-radial domain. This allows the (f> 
dimension to be treated equally to the r — dimensions. The aspect 
ratio (as volume-averag ed within \0 - n/2\ < [ff 1 ],) is given as A, 
in Table [2] where r H is the horizon (focusing on the geometrically 
thinning disk), r, = 20r s (but applies for 50r g > r > 5r e such that 
A is uniform for much of the flow), and r„ = 100r g (showing the 
aspect ratio changes at large radii due to focusing on the jet). The 
MB09D/Q models use r, = 10r s and r = 20r g . 

The MRI is resolved for grid cells per wavelength (Eq. {32}), 

t!.v,MRI = —7 , (43) 

of 2,.MRi > 6, for x = 6,<f>, where A r as dx {l> (dr/dx {i> ), A e as 
rdx (2 \d8/dx^ 2> ), and ss r sin 9dx {3 \d(f>/d/ 3) ). Volume-averaging 
is done as with 5,;,mri, except v Xt \/h x and |Q rot | are separately 6, <f>- 
volume-averaged before forming 2_ v , M ri. Averaged ge.MRi, Q^.mri 
are typically 30% larger than gs, we ak,MRi, 2*,weak,MRi for toroidal 
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Table 2. Numerical Model Parameters 



ModelName 


GridType 


N r 


N 






^out 


A<p 






A Ta 


t=0,(i,o] 


I I 

' f 


AO 94RfN40 


T-TvnF yn 
1 i y L/i_jA ij 


272 


128 


256 


0.855 


26000 


2n 


4.6:1:4.3 


1.1:1.1:1 


1.2:2.1:1 


17 11 


8000-17000 


AO 94RfTM100r1 


T-TvnF yn 
1 i y uijA ij 


136 


64 


128 


0.727 


26000 


2n 


4.6:1:4.2 


1.2:1.1:1 


1.2:2.1:1 


5.6, 3.4 


8000-1 3000 


A0.94BfN100c2 


HypExp 


136 


64 


64 


0.727 


26000 


2n 


4.7:1:8.5 


1:1:1.8 


1:1.8:1.7 


5.8, 3.4 


8000-12000 


AO 94RfNI00c3 


T-TvnF yii 
i j. y uijA lj 


136 


64 


32 


0.727 


26000 


In 


4.6:1:17 


1:1:3.6 


1:1.8:3.4 


5.6 3.4 


8000-1 3000 


A0.94BfN100c4 


HypExp 


136 


64 


16 


0.727 


26000 


In 


4.5:1:33 


1:1:7.2 


1:1.8:6.7 


5.7,3.4 


9000-13000 


A0.94BfN40c5* 


HypExp 


272 


128 


1 


0.855 


26000 


In 


4.7:1:1100 


1:1:230 


1:1.8:220 


16, 11 


8000-13000 


A-0.94BfN40HR 


HypExp 


272 


128 


256 


0.855 


26000 


In 


4.6:1:4.3 


1.1:1.1:1 


1.2:2.1:1 


17, 11 


10000-15000 


A-0.94BfN30 


HypExp 


136 


64 


128 


0.727 


26000 


In 


4.4:1:4 


1.2:1.1:1 


1.2:2.1:1 


11, 6.5 


8000-13000 


A-0 5RfN30 

A U.JDll'l JU 


T-TvnF Yn 
1 1 y uijA ij 


136 


64 


128 


0.735 


26000 


2n 


4.2:1:4 


1.1:1.1:1 


1.1:2.1:1 


9.6 6.7 


8000-1 3000 


AO.OBfNIO 




128 


64 


128 


0.808 


1000 


2n 


2.2:1:2 


1.1:1.1:1 


1.1:2.2:1 


19, 12 


8000-13832 


AO 5R1N30 

AU . .J U 111 J W 


T-TvnF Yn 
i j. y uijA ij 


136 


64 


128 


0.735 


26000 


2n 


4.2:1:4 


1.1:1.1:1 


1.1:2.1:1 


10 7 


8000-1 3000 


AO 94RfN30 

AU . ' i U 11^ JU 


T-TvnF Yn 
1 1 y uijA ij 


136 


64 


128 


0.727 


26000 


2n 


4.6:1:4.2 


1.2:1.1:1 


1.2:2.1:1 


1 1 6.8 


8000-1 3000 


A0.94BfN3()r 


HypExp 


136 


64 


128 


0.727 


26000 


2n 


4.6:1:4.2 


1.2:1.1:1 


1.2:2.1:1 


11,6.8 


8000-13000 


AO 94RnNl 00 


T-TvnF Yn 
i j. y uijA ij 


136 


64 


128 


0.727 


26000 


2n 


4.5:1:4.1 


1.2:1.1:1 


1.2:2.1:1 


7.5 8.4 


1 4000-27048 


A-0.94BtN10 


Exp 


128 


64 


128 


0.797 


1000 


2n 


2.2:1:2 


1.2:1.1:1 


1.2:2.1:1 


< 1 


58000-93944 


A-0.5BtN10 


Exp 


128 


64 


128 


0.806 


1000 


2n 


2:1:1.9 


1.1:1.1:1 


1.1:2.1:1 


< i 


80000-129392 


AO ORtNl 




128 


64 


128 


0.808 


1000 


2n 


2:1:1.9 


1.1:1.2:1 


1.1:2.2:1 




58000-96796 


AO 5RfNl 

AU . .J 1J LI 1 1 VJ 




128 


64 


128 


0.806 


1000 


2n 


2:1:1.9 


1.1:1.1:1 


1.1:2.1:1 


< j 


80000-1 35572 

oyj\jyj\j i _j,j *j / 


AO 94RtNI0 




128 


64 


128 


0.797 


1000 


2n 


2.4:1:2.1 


1.2:1.1:1 


1.2:2.1:1 


< j 


58000-9364R 


AO 94RtNtOHR 

/ V ' f ■ / Til 1 1 ' 1U111\ 


xp 


256 


128 


256 


0.894 


1000 


2n 


2.2:1:2 


1.2:1.1:1 


1.2:2.1:1 


< 1 


58000-801 84 

.J OUUV/ ijyf 1 O^r 


MB09D 


ExpOld 


256 


128 


32 


0.79 


1000 


2n 


1.6:1:12 


1.4:1:11 


1.4:1:10 


4.7, 1.3 


2000-3000 


MB09Q 


FxnOld 


128 


128 


32 


0.816 


40 


2n 


1.6:1:11 


1.5:1:10 


1.4:1:9.7 


< i 


1 500-5688 


A-0.9N100 


TNM1 1 


288 


128 


64 


0.7 


100000 


2n 


2.3:1:7.9 


2:1:6.8 


1.6:1:5.4 


7, 7.2 


11000-20060 


A-0.5N100 


TNM1 1 


288 


128 


32 


0.81 


100000 




1.6:1:8.1 


1.4:1:6.8 


1.1:1:5.3 


7.4, 7.5 


1 1000-16335 


A-0.2N100 


TNM11 


288 


128 


32 


0.81 


100000 


n 


1.6:1:8.2 


1.4:1:6.9 


1.1:1:5.2 


7.6,7.5 


11000-15175 


A0.0N100 


TNM11 


288 


128 


32 


0.81 


100000 


n 


1.6:1:8 


1.4:1:6.8 


1.1:1:5.2 


7.6,7.6 


11000-18585 


A0.1N100 


TNM11 


288 


128 


32 


0.81 


100000 


n 


1.6:1:8 


1.4:1:6.8 


1.1:1:5.2 


7.5,7.4 


11000-20000 


A0.2N100 


TNM11 


288 


128 


32 


0.81 


100000 


n 


1.6:1:8.2 


1.4:1:6.8 


1.1:1:5.2 


7.7,7.6 


11000-19350 


A0.5N100 


TNM11 


288 


128 


32 


0.81 


100000 


n 


1.6:1:8.1 


1.4:1:6.9 


1.1:1:5.2 


7.7,7.6 


11000-19540 


A0.9N25 


TNM11 


288 


128 


32 


0.8 


100000 


n 


1.7:1:8.2 


1.4:1:6.7 


1.2:1:5.3 


16, 15 


11000-17350 


A0.9N50 


TNM11 


288 


128 


32 


0.8 


100000 


n 


1.7:1:8.2 


1.5:1:6.8 


1.2:1:5.3 


11, 11 


11000-14385 


A0.9N100 


TNM11 


288 


128 


64 


0.8 


100000 


2n 


1.7:1:8.2 


1.5:1:6.8 


1.2:1:5.3 


7.9, 7.7 


11000-19895 


A0.9N200 


TNM11 


288 


128 


32 


0.8 


100000 


n 


1.7:1:8.2 


1.5:1:6.8 


1.2:1:5.4 


5.5,5.4 


16000-28620 


A0.99N100 


TNM11 


288 


128 


128 


0.83 


100000 


2n 


1.8:1:4 


2:1:3.3 


1.7:1:2.7 


7.5,7.3 


15000-31400 



field and MB09 models. The t = values (see Table[T} and time- 
averaged values (see Table[8]l are measured at same radii as S^mri- 
Turbulence is resolved for grid cells per correlation length 
(Eq. j3T}), 

CW - (44) 

of Qp.cor > 6, for x = r,6,<p and p = n,l,m, respectively. Oth- 
erwise, modes are numerically damped on a dynamical timescale 
(even Q = 5 would not indicate the mode is marginally resolved, 
because numerical noise can keep Q as 5 at increasing resolution 
until finally the mode is actually resolved - finally leading to an 
increasing Q 3= 6 with increasing resolution ; as seen by Shiokawa 
|et al.|2012) . Reported Q pxoT take 1 /A v as the number of grid cells 
covering the span of A XiCor as centered on: middle of within the 
used radial span for x = r, 6 = n/2 for x = 6 for the "Disk" and 
9 = for x = 8 for the "Jet", and anywhere for x = <p. For A(f> < 2n, 
Q^.mru 2m,cor * is required to avoid truncating the mode (as 
happens in model A0.9N25). 



5 FIDUCIAL THICK DISK MODEL 

Our fiducial model, A0.94BfN40, consists of an initially weakly 
magnetized thick accretion disk around a rotating (a/M = 0.9375) 
BH. The initial magnetic field consists of field loops that alternate 
polarity with radius (the "flipping" field geometry). Model param- 
eters given in Table[T|and Table[2] 

5.1 Initial and Evolved Disk Structure 

Figure [T] and Figure [5] show color plots of b 2 and field line con- 
tours (contours of integrated over (p, so is axially symmetric) 
for the initial and quasi-steady-state evolved solution, respectively. 
The initial solution consists of a radially extended thick torus within 
which several weak field loops (of alternating poloidal polarity) are 
embedded. The disk is geometrically thick with ff 1 ~ 1 and also 
quite thermally thick with cjv mt > 1 and c s /vk ^ 1 through-out 
the solution both initially and at late times. 

The evolved solution, shown in Figure [2] shows that the first 
two field loops have been completely accreted or ejected. During 
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Figure 1. The fiducial model's initial (t = 0) state consists of a weakly 
magnetized geometrically thick torus around a spinning (a/M = 0.9375) 
BH. b 2 is shown as color with legend. Black lines show (integrated over 
all <p). The inner-most-radial field loop is small and shows up as a high 
magnetic energy density with no field lines shown, the next field lines (black 
dashed lines) show the 2nd field loop, and the next field lines (black solid 
lines) show the 3rd field loop. The thickest black solid line is the same field 
line in Figure|2]that tracks the outer-most-angular part of the BH-driven jet. 

accretion of the 3rd field loop, the magnetic flux reaches a maxi- 
mum saturated value over a long period of inflow equilibrium. 

Figure [3] shows an instantaneous meridional snapshot of the 
flow-field, which shows significant circulation. The field lines 
threading the BH are a dipolar split-monopolar magnetosphere. 

5.2 Overall Time Dependence 

Figure [4] shows a typical snapshot for the rest-mass density, field 
lines, and fluxes (M, T, and if) on the BH, through r = 50r g in the 
jet, and at r = 50r g in the magnetized wind. 

The flow consists of a magnetized polar jet and a turbulent 
equatorial disk inflow. Near the BH, the rest-mass density in the 
inflow rises substantially due to vertical compression by the accu- 
mulated polar magnetic flux that surrounds the BH. The large-scale 
polar magnetic flux forms a semi-permeable magnetic barrier to 
the massive inflow, which is forced to undergo non-axisymmetric 
accretion through magnetic Rayleigh-Taylor instabilities. 

The BH's magnetic flux dominates the mass influx with T H ~ 
17 during the quasi-steady-state period. Any additional magnetic 
flux that temporarily accretes onto the hole is ejected in magnetic 
Rayleigh-Taylor modes that push the flux back into the disk. This 
suggests that the magnetic flux near the BH has reached a max- 
imum saturation point via some force balance condition. Because 
T » 1, one expects the BZ effect to be activated, and indeed the en- 
ergy extraction efficiency is high at r) ~ 200%. Most of the energy 
extracted from the BH reaches the jet at large radii (i.e. 77H ~ f/j). 
The shown temporal behavior tracks the fact that 77 oc T 2 . At the 
latest times the efficiency drops as the magnetic flux (from the 3rd 
field loop) begins to be destroyed by an incoming polarity field re- 
versal (outer part of 3rd field loop and inner part of 4th field loop). 



Figure 2. The evolved (; ss 10412r g /c) state of the fiducial model (other- 
wise similar to Figure^ consists of strongly magnetized gas near the BH 
that launches a jet that is visible as the vertical beam of large electromag- 
netic energy density. The jet is nearly cylindrical at small radii due to the 
pressure support from the thick accretion flow. All field lines in the 1st and 
2nd field loops (shown in Figure^ have been accreted or ejected in an out- 
flow. The 3rd (solid black lines) and 4th (dashed black lines) evolved field 
loops are shown, where the thickest solid black line traces the field line that 
connects to the outer angular part of the jet from the BH. 




Figure 3. The evolved (t ss 10412r g /c) state of the fiducial model showing 
velocity flow lines (traces of v-, shown as grayscale) and field lines (A, 
integrated over all cj>, shown as red lines). Also shown are contour lines, 
where — u, = 1 (most particles are unbound in the jet and in the wind at 
larger radii) is a thick black line,/? = 2 (beyond where corona and jet exist) 
is a thick cyan line, and b 2 /po = 1 (beyond where the relativistic jet exists) 
is a thick blue line. The accretion flow is turbulent with circulating eddies. 
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The obvious complete polarity reversal (destruction of inner part of 
2nd field loop by outer 2nd and inner 3rd) occurs at t ~ 2700r g /c 
when M doubles. At f ~ 7000r s /c, magnetic flux is ejected in a 
large magnetic Rayleigh-Taylor disruption allowing M to double. 

Figure[5]shows various quantities vs. time. All quantities are in 
a quasi-steady-state for t > 8000r g /c far after the last field polarity 
inversion at t ~ 2700r g /c and after magnetic flux has accumulated 
near the BH. The mass ejected in the circulating wind (M w „, seen 
as eddies in Fig. [3} dominates the magnetiz ed w ind (M mw . ) and 



jet (Mj) at large radii (r„ = 50r g here); see 5 3.4 for definitions of 
various outflow components. The EM term dominates the MAKE 
term in ?/ H and j H , and the flow has a high efficiency of ^ H ~ 200%. 
The MAKE term is composed of a particle term (i.e. ]f AKE = 1 + 
u t ) and an enthalpy term (i.e. r] EN = u t (u g + p g )/po)- We find that 
^make „ _ 30% as compose( j of rf AKE ~ 63% and t/ en ~ -93% 

(t/ make as -26% in the initial torus, so ?/ EN 67% can be chosen 

for marginally bound inflow). So the inflow is bound as particles 
but fhermo-kinetically unbound due to high enthalpy. 

Accumulation of magnetic flux near the BH leads to geometric 
compression of the inflow as it approaches the horizon. The geo- 
metric thickness drops before the field inversion at t ~ 2700r g /c as 
magnetic flux accumulates and compresses the disk inflow. How- 
ever, during the field inversion, the geometric thickness restores to 
the prior geometric thickness (ff 1 = 0.7) at all radii, which indicates 
that the field (lost during the field annihilation) is responsible for 
the thinning of the dense part of the disk. After the field polarity in- 
version, the magnetic flux re-accumulates near the BH, which leads 
again to the vertical compression of the disk flow. The o--viscosity 
parameter holds steady at about at ~ 0.05. T in the pure inflow 
(;/,. < only) available at large radii (here r = 50r g , giving T outer 
in the plot) is large (the BH and "outer" values are similar for this 
chosen "outer" radius). 

The value of r^ a shows the radius out to which the magnetic 
polarity is the same as on the horizon. As expected, r^ a drops to 
the horizon during the field inversion (destruction of inner part 
of 2nd field loop) at t ~ 2700r g /c. It also gradually drops as the 
next polarity inversion (outer part of 3rd field loop) eats away at 
the magnetic flux outside the BH. The process of field inversion is 
also evident by looking at ThW/Y o (0 (i.e. ratio of time-dependent 
fluxes) corresponding to [the flux on the hole] per unit [flux on 
the hole plus available of the same polarity just beyond the hole]. 
¥n(f)/ x ¥ a (t) ~ 1 is reached during the field polarity inversion, and 
at late times TnW/faW ~ 1 is approached. However, while T 
holds steady, the value of | v P H (0/ l I' a (J)| 1> which indicates that 
much more same-polarity flux is available. This shows that the 
saturated value of T (and so rj) is controlled by some force bal- 
ance condition and not simply limited by initial conditions. Finally, 
WhCO/^h©! ~ 1 shows that the horizon's field is dipolar (/ ss 1). 



5.3 Time-Averaged Poloidal (r - 6) Dependence 

Figure [6] shows the time-averaged flow-field and contours for other 
conditions. The figure is comparable to the snapshot shown in Fig- 
ure [3] The jet region contains significant magnetic flux and same- 
signed polarity field exists near the BH ready to be accreted. In 
the quasi-stationary state, the BH's magnetic flux oscillates around 
its saturated magnitude, whose time-averaged value is determined 
by some force balance condition as managed by non-axisymmetric 
Rayleigh-Taylor instabilities. 
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Figure 5. Quantities vs. time. Top Panel: M out for magnetized wind (M mwo , 
solid line, and mostly middle line), entire wind (M wo , short-dashed line, 
and upper-most line), and jet (Mj, dotted line, and lowest line). Next Panel: 
Efficiency (tjii) for EM (solid line, upper line) and MAKE (short-dashed 
line, lower line) terms. Next Panel: Specific angular momentum Qh) for 
EM (solid line, upper line) and MAKE (short-dashed line, lower line) terms. 
Next Panel: ff 1 at r = {r\\/r g , 5, 20, 100)r ? with, respectively, lines: {solid, 
short-dashed, dotted, long dashed) corresponding to the lowest to upper- 
most lines. Next panel: at at r = i0r g . Next Panel: T on the horizon (solid 
line, mostly upper line) and in the disk at r = 50r g for only the ingoing flow 
(short-dashed line, mostly lower line). Next Panel: ry rl for the radius out to 
where there is the same magnetic polarity as on the hole (solid line). Next 
panel: Magnetic flux on the BH per unit flux available in the flow with the 
same polarity: y ¥u(t)/'Va(t)- Bottom panel: l P tH (f)/ ( l > H(f) ~ 1/', for / mode 
of vector spherical harmonic multipole expansion of A^. In summary, the 
flow has reached a quasi-steady-state at late times. The magnetic flux on 
the horizon has saturated to a large value leading to a high efficiency for 
energy and angular momentum extraction from the BH. 
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Figure 4. Evolved snapshot (see Supporting Information for the movie) of the fiducial model at t » 15612/' g /c showing log of rest-mass density in color (see 
legend on right) in both the z — x plane at v = (top-left panel) and y - x plane at z = (top-right panel). Black lines trace field lines, where thicker black 
lines show where field is lightly mass-loaded. The bottom panel has 3 subpanels. The top subpanel shows M through the BH (Mh), out in the jet (Mj, at 
r = 50r g ), and out in the magnetized wind (M mw , at r = 50r g ) with legend. The middle subpanel shows T for similar conditions. The bottom subpanel shows 
the efficiency (;/) for similar conditions. Horizontal lines of the same colors show the averages over the averaging period, while square/triangle/circle tickers 
are placed at the given time and values. In summary, the efficiency is high at r\ ~ 200%. Also, despite plenty (up to lOx around t ~ 8500r g /c) of same-signed 
polarity magnetic flux surrounding the BH, the magnetic flux reaches a stable saturated value of Th ~ 17 as managed by magnetic Rayleigh-Taylor modes. 
This suggests that the simulation has reached a force balance between the magnetic flux in the disk and the hot heavy inflow. 



5.4 Time-Averaged Radial (r) Dependence 

Figure [7] shows the time-averaged densities, 3-velocities, and co- 
moving 4-fields vs. radius using a density-weighted average to fo- 
cus on heavy disk material. The solution is in inflow equilibrium (3 
inflow times; see section [3/7} only out to r ~ 30r g and has reached 
a single inflow time within r ~ 55-100r g depending upon how one 
defines it. Beyond the BH, the rest-mass and internal energy den- 
sities are quite flat. The rotational velocity is quite sub-Keplerian, 
which is primarily a consequence of thermal pressure playing an 
equal role to total gravity, unlike in thinner disks that are more nat- 
urally Keplerian. This effect is also seen in prior MHD simulations 
Pgumenshchev et al.|200"3j|Pen et al.|2003} . 

The GR viscosity estimate for v r denoted u visc (see above 
Eq. {29}) underestimates the simulation v r when using the a- 
viscosity with total pressure. A more accurate match is found when 
using magnetic pressure. Also, choosing ff* — » ff or ff 1 — > |c s /v rot | 
still leads to a poor fit to v r . Only if we set aiff 1 ) 1 — » 0. 1 at all radii 
does |u v i sc | ~ \v r \ outside the ISCO and inside the inflow equilib- 
rium region. The simulation v r gives e ~ 0.1 for Eq. {T]l. 



Figure [5] shows the fluxes (see section [33} vs. radius as well 
as the field line angular rotation frequency £2 F (using various def- 
initions defined in section |3.6[ l. These quantities are associated 
with conserved quantities such that ratios of total fluxes would be 
constant along flow-field lines in stationary ideal MHD. The total 
fluxes are constant out to large radii, indicating a single inflow time 
is achieved over about 2 decades in radius. The true total fluxes 
are actually even flatter near the BH if one more accurately ac- 
counts for numerical floor injection as done in Tchekhovskoy et al. 
(201 1 1, but this is a small error. Also shown are the components (in- 
flow, jet, magnetized wind, and entire wind) of the mass and energy 
flow. The mass inflow and outflow at large radii follow power-laws 
after sufficient averaging over turbulent eddies. The jet efficiency is 
order 200% and is constant at large radii. 

The winds increase in efficiency with radius, but the jet dom- 
inates the efficiencies of the winds. Power-law fits over the outer- 
radial domain (including the region not actually in inflow equilib- 
rium) for the mass flow rates are M oc r L1 for the inflow and entire 
wind, M oc r° 9 for the jet, and M oc r 0A for the magnetized wind. 
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Figure 6. Flow-field as in Figure [5] except time-averaged and zoomed-in 
on the BH region that is in inflow equilibrium (3 inflow times). The colored 
(green, black, cyan, and blue) thick lines correspond to time-averages of 
quantities Q = {u r , w,,p\ £> 2 /Pol, respectively, at values [Q], = {0, — 1,2, 1}. 
While near the BH the flow has [b 2 /po], > 1 as averaged directly, the dense 
inflow has fo 2 /po < 1 at all radii. The inflow occurs in geometrically thin 
streams not accounted for when computing the unweighted average. In sum- 
mary, the BH is threaded by ordered magnetic flux, and the flow exhibits 
equatorial asymmetry over many inflow times. 



g , the entire wind component has an average velocity of 



At r ~ 30 ; 
order v w /c ~ 0.01 or v w ~ 3000km/s. 

The specific magnetic flux T measures the total absolute ra- 
dial flux. The figure shows that the total radial flux is much larger 
than the inflow-only component because the jet harbors most of the 
magnetic flux. The field line angular frequency fl F ~ fl H /4 (as in 
BZ77's paraboloidal model) in the disk+corona+wind (i.e. "fdc" 
averaging, for full flow except the highly-magnetized jet). 

Figure [9] shows the time-averages for the disk's geometric 
half-angular thickness (ff 1 ), the thermal half-thickness (ff, using the 
density-weighted average), flow interface angular locations, reso- 
lution of the MRI wavelength, approximate a viscosity parameter, 
and magnetic fluxes vs. radius. The magnetic field compresses the 
disk leading to decreasing ff 1 with decreasing radius. While ff 1 <K 1 
near the BH, ff = arctan (c s /v rot ) > 1 and arctan (cJv K ) > 1. 
Hence, the flow is not in vertical hydrostatic equilibrium due to 
the strong magnetic field. Note that using u K instead of v rat in the 
expression for ff gives a similarly large thermal half-thickness. The 
disk-corona and corona-jet interfaces trace the path of the well- 
collimated jet out to large radii. The <2e,MRi » 6 as required to 
resolve the MRI ( Sano et al. 2004), while in the inflow equilibrium 
region S^mri ^ 1/2, indicating the MRI is suppressed. Even using 
SI — > £1 K in SVmri leads to r SdMR1=l/2 « 8r s , so the MRI suppres- 
sion near the BH is not only due to sub-Keplerian motion but also 
disk compression by the accumulated magnetic flux. The horizon's 
time-averaged radial absolute magnetic flux is O r H ~ 200 and there 
is an additional T cq ~ 100 same polarity magnetic flux available. 



5.4.1 Comparison with Gammie (1999) Model 

Figure [TO] shows a comparison between the fiducial model and the 
Gammie ( 19991 model of a magnetized unstratified accretion flow 




0.01 i 



100 



Figure 7. The time-angle-averaged densities, 3-velocities, and 4-field 
strengths using a density-weighted average. Top panel shows rest-mass den- 
sity (po) as black solid line, internal energy density (u g ) as black short- 
dashed line, and magnetic energy density («/,) as black dotted line. Middle 
panel shows negative radial velocity (— v r ) as black solid line, rotational ve- 
locity (v rot ) as black short-dashed line, Keplerian rotational velocity (hk) as 
blue short-dashed line, and a-viscosity theory radial velocity (t! v ise) when 
using pb in denominator for a = a\, (blue solid line) and when choosing 
a fixed ciiff 1 ) 1 =0.1 (blue dotted line). Bottom panel shows comoving 4- 
field spatial components with r, 8, and ip components shows as solid, short- 
dashed, and dotted black lines, respectively. The vertical red line marks the 
ISCO. Vertical solid cyan lines show range from r = I2r g to 3 inflow times. 
The short-dashed vertical cyan line marks a single inflow time. In summary, 
dense filaments maintain Ub/po i 1 and /J » 1 near the horizon, while a 
similar plot (see bottom-left panel in Figure |10) averaged over the entire 
disk+corona+ winds shows ut/po » 1 and fl <s 1 near the horizon. Also, 
the rotational velocity is quite sub-Keplerian. 



within the ISCO. A value of T as 8.7 was used to match the Gam- 
mie model's value of the EM component of the specific angular 
momentum at the horizon, while the effective disk thickness was 
varied to best fit the magnitude of u b at the horizon. All quantities 
use "fdc" type averaging focusing on all parts of the flow except the 
highly-magnetized jet. Compared to the densities vs. radius shown 
in Figure [7] here clearly u b » p on the horizon due to averaging 
over the entire disk+corona+winds. The densities monotonically 
increase towards the BH as the inflow is vertically compressed by 
magnetic field. 

If T is only integrated across the time-averaged disk where M 
is non-zero (i.e. "fdc" averaging: full flow except the highly mag- 
netized part of the jet), then we find T « 1.5. However, within 
the heavy dense filaments (i.e. density-weighted average), we find 
Y ss 0. 1 . Also, j} » 1 in the dense flow, while ft <k 1 over the 
disk+corona+winds. This shows that the heavy filamentary parts 
of the disk inflow are quite under-magnetized compared to their 
surroundings, which is as expected for accretion through efficient 
magnetic Rayleigh-Taylor instabilities. 
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Figure 8. The time-averaged angle-integrated fluxes. From top to bottom, 
panels are: Total mass accretion rate (M), inflow rate (M m ), jet outflow rate 
(Ms), magnetized wind outflow rate (M mw ), entire wind outflow rate (M w ), 
total specific energy accretion rate (£/Mh), efficiency for the jet (solid line) 
magnetized wind (short-dashed line) and wind (dotted line), total specific 
angular momentum accretion rate (j = j/Mu), specific magnetic flux T for 
the total flow (solid line) and pure inflow with u r < (short-dashed line), 
and field line angular rotation frequency per unit BH angular frequency 
(Af/Ah) f° r time-averaged versions of Ap (solid line), A p (short-dashed 
line), Ap (dotted line), |Ap| (long-dashed line), and |flp| (dot- short-dashed 
line). These flp are averaged within the disk+corona part of the flow. Power- 
law fits for mass inflow and outflow rates are shown as short-dashed lines. 
In summary, inflow equilibrium is achieved over a couple decades in radius, 
and the mass outflows follow a power-law behavior. 
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Figure 9. Other time-angle-averaged quantities. From top to bottom, panels 
are: density half-angular thickness (ff 1 , solid line) and thermal half-angular 
thickness (0' , short-dashed line), disk-corona interface angle (8 dc , solid line) 
and corona-jet interface angle (6 CJ , short-dashed line), number of cells per 
fastest growing MRI wavelength (Qa.mri. solid line ; Qe, w eak,MRl. short- 
dashed line) with Qe.MRl = 6 shown as dotted line above where the vertical 
(8) MRI is resolved, number of cells per fastest growing MRI wavelength 
(2<6,MRI, solid line ; Q#,weak,MRl> short-dashed line) with Q^mri = 6 shown 
as dotted line above where the azimuthal (<p) MRI is resolved, number of 
fastest growing MRI wavelengths across the full disk thickness (S^mri, 
solid line ; Srf, w eak,MRl> short-dashed line) with S^mri = 1/2 shown as 
dotted line below where the MRI is suppressed, viscosity parameter (a/,, 
solid line ; atb^g, short-dashed line), radial absolute magnetic flux (<D), and 
equatorial magnetic flux (Tcq). In summary, the disk is compressed by the 
horizon's magnetic flux leading to smaller 3 d as r drops despite little change 
in 8'. Also, the linear MRI is suppressed even in the dense inflow. 



5.5 Time-Averaged Angular (9) Dependence 

Figure[TT]is similar to Figure[7]but for quantities vs. 6 at four dif- 
ferent radii. This also highlights how the disk flow is compressed 
as it approaches the horizon. The time-averaged density is well-fit 
by a Gaussian with width ff 1 ss 0.12 at r = ru even though the time- 
average of the instantaneous ff 1 ~ 0.06. This is because the time- 
average of density blurs the width of the narrow disk that oscillates 
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Figure 10. Comparison between the accretion flow for the fiducial simula- 
tion ("fdc" averaging, focusing on entire flow except the highly-magnetized 
jet, shown by solid lines) and the Gammie 1 1 999 1 model for an unstratified 
magnetized inflow within the ISCO (shown by dashed lines). In all panels, 
red vertical lines shows the location of the ISCO. Top-left panel: Shows 
the radial 4-velocity, where the Gammie solution assumes «' = at the 
ISCO. Finite thermal effects lead to non-zero i/ at the ISCO for the simu- 
lated disk. Bottom-left panel: Shows the rest-mass density (po, black line), 
the internal energy density (u g , magenta line), and magnetic energy density 
(iq,, green line). Top-right and bottom-right panels: Show simulation results 
and Gammie solution for the specific particle (PA) flux (black line), specific 
enthalpy (EN) flux (magenta line), and specific electromagnetic (EM) flux 
(green line). Horizontal dotted black lines are the Novikov-Thorne thin disk 
solution (i.e. particle/radiation term). Thi s figure is comparable to fig. 10 for 
a somewhat thick (¥ ~ 0.2-0.25) disk in McKin ney & Gammie|(2004} and 



to fig. 1 1 for a thin ( 

thick disk simulations, the model fits j and b l but not E. 



0.05) disk in Penna et al. 

i7 



(2010) . As for the prior 




in height about equal to its own height. Since the time-averaged 
density is a Gaussian, the estimate of the thermal thickness p given 
by Eq. {Bj would have given 0' p ~ ¥ if the disk were in hydrostatic 
equilibrium. That ff' <k 6' close to the BH (as shown in Figure [9} 
shows that the thermal content remains relatively unchanged de- 
spite significant geometric compression by the polar magnetic flux. 
The figure also shows that the total pressure is roughly constant 
with angle as dominated by the magnetic pressure. The behavior of 
v e , v$ near the polar axes is affected by the numerical floor mass 
injection, although this region contains little energy or energy flux 
compared to other angles. 

Figure [l 2| shows the horizon's value s of quantities related to 
the BZ effect (Blandford & Znajek 19771. The simulation's fluxes 
are computed via Eq. |l0|. The "full BZ-type EM formula" referred 
to in the figure uses the EM energy flux computed from equation 33 
in McKinney & Gammie (2004), which only assumes stationarity 
and axisymmetry (rather than also small spin in BZ77) and uses the 
simulation's Slp(6) and B r (9) on the horizon. This figure shows that 
most of the horizon is highly magnetized due to accretion occurring 
through a magnetically compressed inflow. 

The agreement between the simulations and the BZ picture is 
excellent for the highly magnetized regions, where roughly £ip ~ 
fin/4 near the disk-jet interface (here, fl F is the time-average of 
Eq. f27}). While the simulation is roughly consistent with BZ's 



Figure 11. Similar quantities as in Figure [7] except plotted vs. 8 at 
f — {rii/>'g, 4, 8, 30}r g (respectively: solid, short-dashed, dotted, and long 
dashed lines). If numerical density floors were activated at some space-time 
point, then po = u g = was set there. This shows up even in these time- 
averaged densities as a sharp drop-off, because the jet maintains high mag- 
netizations at the poles during the interval time-averaged. In summary, the 
disk is compressed by the polar magnetic flux as the inflow approaches the 
horizon, and the total pressure is roughly constant with 8. 

paraboloidal solution, the equatorial fi F is somewhat suppressed 
due to the disk inflow. Also, near the polar axes, Op is affected by 
ideal MHD effects and numerical floor mass injection. In the last 
~ 3 grid cells, limited resolution affects Sip. However, both the grid 
and jet collimate, which leads to a well-resolved £1? near the axes at 
slightly larger radii. Also, the most energetic part of the jet is near 
the disk-jet interface that dominates the dynamics at all radii. Even 
if we used a paraboloidal extension of our simulation data into the 
region where the numerical floor injection is important, this only 
changes the total efficiency by less than 20%. 

5.6 Azimuthal (0) Time Dependence 

Figure [T3] shows a large magnetic Rayleigh-Taylor disruption at 
t « 15172r g /c once T has reached its saturated value and the 
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Figure 12. Time-0-averaged quantities and flux integrals on the horizon 
as a function of 8. From top to bottom: 1) Field line rotational angular fre- 
quency (£Jf/£Jh) for simulation (solid line), 1 st-order-in-spin accurate value 
for monopolar (short-dashed line) and paraboloidal (dot-long-dashed line) 
BZ solutions ; 2) Rest-mass flux (Mu) ; 3) Electromagnetic (EM, solid line) 
and matter (MA, dot- short-dashed line) efficiency (i7h), along with the full 
BZ-type EM formula without any renormalization (dotted line) ; 4) Elec- 
tromagnetic (EM, solid line) and matter (MA, dot- short-dashed line) spe- 
cific angular momentum flux (ju = j^/Mu), along with the full BZ-type 
EM formula without any renormalization (dotted line); 5) Gammie param- 
eter (Th) for the simulation (solid line), and the BZ model for the cases: 
Oth-order-in-spin accurate monopolar field (short-dashed line), Oth-order- 
in-spin accurate paraboloidal field (dot-long-dashed line), and 2nd-order- 
in-spin accurate monopolar field (long-dashed line). These BZ versions are 
normalized so total magnetic flux is the same as in the simulation. Notice 
how the 2nd-order-in-spin accurate monopolar BZ model fits the simula- 
tion result quite well. For the last 3 panels, the divisor is (implicitly) Mh 
that has been fully angle-integrated to a single value. So, rjw, ju, and Th 
show the angular dependence of .Eh, Jh, and "Fh, respectively. In summary, 
the agreement between the simulation and the BZ picture is excellent. 



flow can no longer accept new magnetic flux on the BH. While an 
atypical looking outburst, other more typical snapshots, such as at 
t « l\212r g lc, show roughly similar behavior including the same 
range of density variations. Such outbursts occur because magnetic 
flux is temporarily added to the BH and exceeds the saturated value, 
but then that extra flux is ejected back in magnetic interchange 
modes. Here, \m\ = 1,2 modes appear to dominate. 

Figure [14] shows the density vs. <p at three radii (r = 
ra, 4r g , 8r s ) for t ■ 
the average with weight w 



100 r 



15172r g /c (as shown in Figure 13 I. p shows 
1, and p^ q shows the density for a 
single grid cell cut exactly at the equator. Typical density variations 
averaged over the flow are up to factors of 10 near the horizon, 
while at any given 8-<p angles (here the equator) the density varies 
considerably (down to the limits of the numerical density floors). 
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Fi gure 14. Shows rest-mass density at t ss 15172r„/c (as shown in Fig- 
ure 131, as 0-averaged over the entire flow (p'q S ) and as cut along the equa- 
tor (p^ 1 ) at three radii (in all panels: r = ru shown as solid line, r = 4r g 
shown as short-dashed line, and r = 8r„ shown as dotted line). The drop- 
outs in density correspond to where the mass is dominated by numerical 
floor mass injection due to the high ut/po in such regions, so we have set 
po = there to show where the density could be smaller than the floor. In 
summary, the variations in density become quite large near the BH. 



This shows how the dense part of the flow is forced through the 
magnetosphere via magnetic Rayleigh-Taylor modes. 



5.7 Time-Averaged Azimuthal (0) Dependence 

Figure [T5] shows the normalized time-average of the Fourier de- 
composition in the (^-direction ([|a,„|]//[|ao|]r) using Eq. fiO\ for 



the "Disk" and "Jet" with quantity Q as pou r , po, u g , b 2 , and T EMr t 
that gives, respectively, the actual M, comoving mass M , comov- 
ing thermal energy E g , comoving electromagnetic energy E B , and 
actual electromagnetic energy flux E EM . As with flux ratios, we 
compute []a m ]]r/[|ao|], instead of [|a m /aolL because the latter would 
exaggerate mode power during transient moments when «o(0 be- 
comes small and potentially zero. 

We find that there is large \m\ > power in M and £em, but 
there is little power in |m| > for b 2 in the jet. Since |a m >o| ~ \a m= o\ 
for M, this shows that accretion occurs primarily through non- 
axisymmetric modes. Also, |a m> ol ~ |« m =ol f° r the jet electromag- 
netic power even at r = 8r g , so the jet power contains significant 
non-axisymmetric structure. The \m\ = 1 dominates all |m| > 0, but 
\a m \ and |ai| are similar up to m ss 20. 

Figure [15] also shows the azimuthal correlation length scale 
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Figure 13. Snapshot at t ss 15172r g /c for the fiducial model during one of the larger magnetic Rayleigh-Taylor disruptions, just after which the mass accretion 
rate increases. Otherwise similar to the upper panels in Figure[4] During this disruption, \m\ = 1,2 modes appear to dominant the non-axisymmetric accretion. 



(Eq. {ST}), which is always well-resolved beyond the horizon. Note 
that if we compute Eq. ^30\ as a spectrum of the averaged flow 
rather than as an average of the spectrum, then even on the hori- 
zon the correlation lengths are very well-resolved. The correla- 
tion length scale A,/, XOI is related to the mode m cor = 2njX^ COI . 
For example, for M (i.e. the density in the disk), we find m cor « 
33.8,15.9,13.2,8.5 at r = r H ,4r g , 8r g , 3f> g , respectively. With 
¥ ~ 0.06,0.13,0.29,0.59, one obtains A^Jff 1 » 3.1,3.0, 1.6, 1.3, 
respectively. So as the magnetic flux near the BH vertically com- 
presses the inflow, the azimuthal correlation length drops, but not 
in proportion to the disk height. 



5.8 Quasi-Periodic Oscillations 

BH accretion disks are observed to have QPOs roughly classified 
as either high-frequency QPOs (HFQPOs) or low-frequency QPOs 
(LFQPOs). In BH x-ray binaries, QPOs are seen in only specific 
states (Fender et ah[|2004l |Abramowiczl[2005] |Remillard & Mc-] 
|Clintock||2006[ >, such as the steep-power law (SPL) state (related 
to the very high state and intermediate state, during which power- 
ful transient jets are observe d) ([Done & Gierliriski 2003| |Gierlinski| 



& Newton 2006; |Oda et al.||2009| >. For BH x-ray binaries, HFQ 
POs range from 100Hz (GRS1915+105) to 300Hz (1655-40) cor- 
responding to a period of r ~ 70-130r g /c. HFQPO frequencies are 
sometimes in a 3:2 ratio ( TPsaltis et al.|l9 99 ). QPOs are not expected 
to be as easily seen in AGN as for x-ray binaries ( Vaughan & Utt 



ley|2005|l; but some AGN may have QPOs, such as SgrA* l Meyer 



et al. ||2008 1 |Dolence et aTJ|20T2l and RE J1034+396 ([Gierlinski 
et al.|2008| >. NS QPOs have similar features (van der KIis|1998> . 

BH x-ray binary QPOs are not seen (or are very weak) in the 
high-soft state that is typically associated with a thin disk domi- 
nated by MHD- MRI turbulence (|Wijnands & van der Klis|1999| 
|Remillard||2005l |Rubio-Herrera & Lee||20051 |Machida & Mat-| 
|sumoto||2008l |Mondal et aTJ|2"009"l >. The lack of QPOs in such a 
state is expected because non-linear turbulence tends to reduce co- 
herence. Also, the existence of a bright hard state defies explana- 
tion by standard viscous models that would predict the bright state 



should be soft ( |Oda et al.|20 09 ). This indicates that a qualitatively 
new accretion state (like the MCAF state) may be required. 

Analytical models and simulations have demonstrated QPOs 
through various mechanisms. The hope is that QPOs could be used 
to measure BH spin and mass (Stella & Vietri 1999). Disk mode 
oscillations may cause either HF or LF QPOs |Reynolds"& Miller 
[2009l|O'Neill et al.|2009|[20TT) , disk precession causes LFQPOs in 
GRMHD simulations of tilted disks pexter & Fragile|20lT) , and 
dynamo field oscillations in local MHD simulations show LFQPOs 
( [Davis et al.|20"T0l |Guan & Gammie|20lT) . In some cases QPOs 
may be seen, but their statistical significance is highly uncertain. 
This includes HFQPOs seen in pseudo-Newtonian MHD simula- 
tions ( |Kato|2004) . GRMHD simulations have shown HFQPOs in 
tilted disks (Henisey et al. 2009). However, in the latter case higher 
resolutions were found to eliminate the QPOs as distinct features 
because MRI turbulence destroys their coherence (Henisey, 2010, 
priv. comm.). Other GRMHD simulations with radiative transfer 
show some HFQPO features ( Schnittman et al. 2006 1, but their sig- 
nificance is highly dependent upon model of "continuum" spectrum 
as a single power-law, while their power spectrum is identified as a 
broken power-law. So extra apparent power appears near the break. 

Our GRMHD simulations also show coherent HFQPOs. We 
identify the coherence of the QPO by its quality factor Q « v/(Av) 
for frequency y, where Ay is the full width at half of the maxi- 
mum power (FWHM) (e.g. if the QPO distribution is modelled 
by a Lorenzian). All our thick disk poloidal field simulations at 
both resolutions (136 x 64 x 128 and 272 x 128 x 256) show co- 
herent QPOs during some part of the simulation at various radii 
(r = 7'h, 4r g , 8r g , 30r g ), while the QPOs are most coherent for the 
high-resolution fiducial model. This expected resolution depen- 
dence on Q implies that damping at low resolutions can make it 
difficult to resolve a coherent QPO. Radiative transfer ( |Broderick| 
|& McKinney|2010|[Shcherbakov et al.|2010| [Dexter et al.|2011| > 
is necessary to see if these QPOs are observable. In this paper, we 
only consider the MHD dynamical properties of the HFQPO. 

Figure [T6] shows b^, averaged over (f> = to = n/4 (with no 
weighting) vs. 9 and time t at r = 4r g (averaged over ±0.4r g with 
no weighting) for the fiducial model. This shows possible QPOs. 
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Figure 16. vs. t and 8 at r = 4r g . The field polarity switches at 
t ~ 2700r e /c, after which magnetic flux accumulates and saturates. Af- 
ter saturation, the jet-disk QPO (JD-QPO) becomes coherent, here visible 
as horizontal stripes at late times. 



Figure 15. Time-averaged Fourier amplitude for m modes. Each panel 
shows result for r = ru, 4r g , Hr g , 30r s (respectively, solid, short-dashed, 
dotted, long-dashed lines), except the top panel that only shows results 
for r = r\\,4r g with same line types. The power in all |m| > modes 
relative to the \m\ = mode is shown for each panel respectively for 
radii r = ru,4r g ,8r g ,30r g . For each line, the vertical bar corresponds to 
»26cells = where m > m^ ce \\ s are numerically damped on a dynamical 

time. For each line, squares mark the correlation length's m = m cor mode, 
so azimuthal structure (e.g. turbulence) is resolved if m cor < m6 Ce ils. In sum- 
mary, azimuthal structures are well-resolved across all quantities within the 
causally-connected region outside the horizon. Changes in m col with radius 
in the "Disk" partially track changes in the disk thickness. Most power is 
at \m\ ~ 1, but significant power extends up to |m| ~ 20 before dropping 
off more rapidly. Because |a,„>o| ~ \a m -a\ for M, mass accretion is highly 
non-axisymmetric (see also Figure [14) . 



Figure[T7]shows a spectrogram for the power density vs. time 
and frequency for b 2 at r = 4r g at 9 ~ 2.9 (i.e. deep within the 
jet). Power density is shown in commonly-used units: an integral 
of power density over a frequency range gives back the square of 
the fractional root-mean-squared (rms) amplitude of the variability 
in the original time series. 



Figure 18 shows the Fourier transform of b-(t) at the equato- 
rial plane (at 9 = 0), in the disk (at 6 = 1 [9 d ]„ where [9*], is the 
time-averaged 9 d ) and in the jet (at 9 ss 2.9) at r = 4r g as averaged 
from = 0-7T/4 to reduce noise while allowing us to resolve sev- 
eral \m\ modes. The power density is given in units of (rms/mean) 2 
multiplied by the total period (t = 12000-16000r g /c, focusing on 
the high-Q period) over which the Fourier transform is performed, 
as a function of frequency (in units of c/r g ). We find that Q ~ 100 
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Figure 17. Spectrogram for the power density vs. time and frequency for 
b 2 (with legend, where white is below lowest power in legend) at r = Ar g 
at 9 w 2.9 (i.e. deep within the jet). This shows the high Q oscillation starts 
once flux saturates for t > 5000r s /c, where ~ 3 jet-disk QPO (JD-QPO) 
modes are visible. The frequency changes little as M varies. 



in the jet, Q ~ 10 one scale height above the disk plane, and Q < 10 
in the disk plane. Other radii (e.g. r - r H , 8r s , and 30r s ) in the jet 
also show similar QPOs as features move out in the jet. 

What is the nature of the HFQPOs in our simulations? Once 
poloidal magnetic flux has accumulated and reached its natural sat- 
urated limit, a semi-permeable magnetospheric barrier forms be- 
tween the heavy disk inflow and magnetic flux threading the rotat- 
ing BH and its jet. This magnetospheric interface is where mag- 
netic field and density change somewhat abruptly. Linear stabil- 
ity analyses of such magnetospheric interfaces predict them to be 
unstable to Rayleigh-Taylor and Kelvin Helmholtz modes, which 
drive QPOs at spherical harmonic m modes based upon some ro- 
tational frequency (SI) such that Qqpo ~ mSl \hi & Narayan|2004) 
Fu & Lai 20121. These linear stability analyses predict that low m 
modes dominate the non-linear dynamics and that QPOs should ap- 
pear. However, these models cannot determine what controls SI, the 
power of the m modes, or the coherence of the QPO. 

Our MCAF simulations validate the prior linear stability anal- 
ysis. We confirm that the magnetospheric interface is indeed unsta- 
ble and that low-|m| (primarily \m\ = 1) modes dominate. We also 
show that the unstable jet-disk interface drives coherent HFQPOs. 
A movie of snapshots like Figure [4] shows that once poloidal mag- 
netic flux saturates, the jet and disk oscillate together. 

Why do such QPOs appear in MCAFs and not MRI- 
dominated accretion flows? QPOs seen in prior MRI-dominated 
accretion flow simulations show (e.g.) m = 1 spiral disk modes ex- 
tending out in radius (Dole nce et al.|2012] >. Other QPOs ( |Henisey| 
|et al.||2009) were absent at high resolutions due to decoherence 
by properly-resolved MRI-driven MHD turbulence (Henisey, 2010, 
priv. comm.). Such MRI-dominated simulations did not reach sat- 
uration of poloidal magnetic flux because of their limited initial 
poloidal flux. Yet, only once the poloidal flux has saturated does 
the magnetospheric interface form. In our models, the QPOs are 
driven due to the presence of an unstable interface between the jet 




10 10 

Frequency [/, c/r ff l 

Figure 18. Power Density for b 1 at r = 4r g in the disk plane (top panel), one 
geometric half-angular thickness ((/') above the disk plane (middle panel), 
and deep within the jet at 8 a 2.9 (bottom panel) shown as black lines. A 
smoothed (over 10 frequencies) version is shown as red lines. The disk has 
QPOs with quality factor Q < 10, a disk scale-height above the disk plane 
has Q ~ 10, while the jet itself shows a few different harmonics with up to 
Q ~ 100. The jet harmonics correspond to \m\ = {1,2, 3} modes based upon 
the field line angular frequency (driven by the black hole rotation frequency) 
at the disk-jet magnetospheric interface where the jet-disk QPO (JD-QPO) 
mechanism operates. 
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and disk where force balance has been achieved between magnetic 
forces by the jet and (e.g.) ram forces by the disk inflow. 

What sets the frequency of the QPOs in our simulations? In 
prior magnetospheric QPO models, the QPO is driven as an inter- 
face instability that oscillates at an m-mode-based rotational fre- 
quency £1. The magnetospheric interface in our simulations is be- 
tween the jet and the disk. The strong field at the disk-jet inter- 
face forces the plasma rotation frequency (Q.) to be similar to the 
field line rotational frequency £1 F , so all Q are similar at the in- 
terface. The QPO frequency is thus set by the BH spin frequency 
(fi H - a /(2Mr H )) due to the BH dragging the field at angular fre- 
quency Qp as Q.n/4 at the jet-disk interface, as consistent with BZ's 
paraboloidal solution (see Figure [T2|. 

The QPO frequencies seen in Figures |17|18| correspond to 
|m| = 1, 2, 3 modes at the jet-disk interface. That is, flopo * mO.? « 
mQu/4. The \m\ = 1 mode dominates, as shown in Figure [T5] So, 
for example, the dominant QPO is due to the |m| = 1 mode that in 
the fiducial simulation has period r ~ 70r g /c, which agrees with 
t » lit /(mQ.^/4) ~ 123r g /c for m = 1 and a = 0.9375 from 
the above analysis. Other BH spins (e.g. \a/M\ < 0.4) might satu- 
rate at Qqpo ~ mil at the (e.g.) ISCO due to the disk dominating 
the plasma rotation rate near the disk-jet interface for such models 
( McKi nney & Gammie|2004| |McKinney|2005fr . 

In summary, the period r ~ !Qr g /c for a/M ~ 0.9 for our 
\m\ = 1 mode (with longer periods expected for lower BH spins) is 
consistent with the range of HFQPOs observed in BH x-ray bina- 
ries. Also, the disk-jet interface (driving the QPO) harbors a large 
electromagnetic energy density, so the disk-jet interface can domi- 
nate (non-thermal) synchrotron emission. Our work shows that co- 
herent HFQPOs may be initiated by large-scale magnetic flux near 
the BH, but more work is required to test their observability. 



6 DEPENDENCE UPON FIELD GEOMETRY AND 
STRENGTH, SPIN, AND RESOLUTION 

This section provides results for all our models. We also take this 
opportunity to tabulate results from the brief letters by |McKinney &| 
|Blandford|p009) ; |Tchekhovskoy e t al. (2011}. Quantities are time- 
averaged from r? to T" given in Table |2| The chosen T" is at least 
after the dimensionless fluxes through the horizon have reached a 
quasi-steady-state. Typically T" = Tf. Units in tables are GM = 
c = r g = 1. Such tables are computed in other works ((McKinney &| 
|Garnmie|2004||De Villiers et al.|2005||Hawley & Krolik|2006| >, and 
our tables should prove useful for comparisons with future works. 

Summarizing our MB09 models: Results are similar to most 
prior MHD simulations of moderately thick disks. The flows are 
dominated by the local MRI, nearly Keplerian, and not efficient. 
For MB09D, T". = 3000, after which turbulence decays. The model 
MB09D focused on the large-scale jet that is only moderately af- 
fected by how the disk behaves at late time, and the jet's correla- 
tion lengths are well-resolved beyond the horizon. Unlike our other 
models, a color plot (not shown) of b<p or b 2 vs. t and 6 shows a but- 
terfly type LFQPO behavior with a period of roughly 15-30 times 
the orbital period at r = 4r g or 8r s , as similar to fig. 10 in |Guan &| 
|Gammie| ( |2%TT] and figs. 8 and 12 in |Davis et aLlpOlOt . 

Summarizing our "thinner disk" TNM11 models: Results are 
similar to the thick disk poloidal field models. Plotting (not shown) 
all prior figures, some similarities and differences are notable. The 
TNM1 1 models are more Keplerian due to ff < 1. The efficiencies 
are 7, ~ 100% for a/M > 0.9 and ~ 30% by r = 100r s 

for a/M = 0.99. On the horizon, Q. F /S1 H vs. 6 more closely fol- 



lows the profile expected for a paraboloidal field due to the lower 
numerical density floors. S^mri 5; 1/2 inside r ~ 20r s (dependent 
upon initial and simulation duration), indicating that the MRI is 
suppressed over smaller radii. The value of a at r ~ 10r g is larger. 
Compared to the thick disk models, the density-weighted inflow is 
even more weakly magnetized (i.e. b 2 /po <k 1 and /? » 1). Also, 
e ~ 0.05-0.1, as used in Eq. {TJ. No coherent HFQPOs are seen in 
the thinner disk models, but the quantities (e.g. rjn) are much more 
variable. The lack of HFQPOs and the higher variability may be 
due to the disk being thinner. However, for some of these models, 
this could also be due to lower Ng, N4, available for capturing mag- 
netic Rayleigh-Taylor modes, e.g., the azimuthal correlation length 
is only marginally resolved with 2„,_ col < 6. (However, the time- 
averaged quantities are actually well-resolved.) 



6.1 Disk Thickness, Disk-Corona, and Corona- Jet Interfaces 

Table pi shows the number of grid cells (Ng H ) at the horizon that 
cover the angular span of the disk's geometric half-angular thick- 
ness (ff^) (ff 1 is computed by Eq. (|12[), the disk thickness at the 
horizon (ff^) and other radii r (dfj r ), the thermal half-angular thick- 
ness (ff 20 is ff at r = 20r g as computed by Eq. d 1 3b), the disk-corona 
interface angular locations at the same radii (pf^, defined by where 
P = 1), and the corona-jet interface angular locati ons a t the same 
radii (£& , defined by where b 2 /p Q =1). See section 3.4 for details. 
In summary, jets in the thicker disk models collimatebetter. 

The rest-mass density geometric half-angular thickness (ff 1 ) 
is quite different than the thermal half-angular thickness (ff = 
arctan(c J /w ro t)), the latter being an estimate of ff 1 if gas pressure 
balances vertical gravity. Close to the BH, ff 1 <k ff is controlled 
by magnetic forces, and ff ~ (r only at large radii. About 6 cells 
should span the full disk to resolve it at a basic level, and most 3D 
models satisfy this even on the horizon (~ 20 cells for the fiducial 
model). The "magnetospheric radius" (r m ) is not a sharp boundary, 
and instead magnetic forces gradually compress the disk. 



6.2 Mass Accretion and Ejection Rates 

Table[4]shows M computed via the first row of Eq. 1 14 1 through the 
BH horizon (M H ), inflow-only mass accretion rate (M m ) at an inner 
radius of r, = \0r g giving Mj„.j and an outer radius of r„ = 50r g 
(except model MB09Q that uses r = 30r g due to its limited radial 
domain) giving M m o . Also shown are M through the jet (Mj), the 
magnetized wind at r, (M mw> ;) and r D (M mw , a ), and the entire wind 
at the same radii (M w j and M Wi0 )- 

Most models show massive winds, and the poloidal models 
show some mass carried in the magnetized wind and jet. Our mod- 
els are in inflow-outflow equilibrium with M- m „ - M H ~ M w o + Mj 
(i.e. all outflowing material, that is not in the jet, is part of the "en- 
tire wind"). Little mass is accreted in the 2D poloidal models be- 
cause magnetic flux accumulates and fully (except through periodic 
reconnection events) suspends the inflow. 



6.3 Energy Efficiency of Hole, Jet, and Winds 

Tab le [51 shows the flow efficiency computed via Eq. \17\ (see sec- 
tionpM), where 77 = rf M + t/ make and 77 MAKE = i] mKE + rj m . These 
efficiencies are computed at the horizon (^h)> f° r the jet (rjj), for the 
magnetized wind (?/ mw ), and for the entire wind (77 w ). The winds 
are measured at r, and r„ given above, while the jet is measured 
at r . The radial asymptotic efficiency is r] m ~ rjj + r] mvlo , because 
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Table 3. Grid Cells across Half- Thickness at Horizon, Half-Thicknesses of Disk, and Location for Interfaces for Disk-Corona and Corona- Jet 
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0.13 


0.3 


0.51 


0.47 


0.053 


0.097 


0.31 


0.45 


0.072 


0.22 


0.88 


1.3 


A0.2N100 


4.2 


0.057 


0.12 


0.29 


0.49 


0.51 


0.053 


0.1 


0.3 


0.25 


0.068 


0.21 


0.82 


1.3 


A0.5N100 


3 


0.042 


0.11 


0.26 


0.51 


0.46 


0.041 


0.1 


0.32 


0.36 


0.055 


0.21 


0.79 


1.3 


A0.9N25 


5.2 


0.068 


0.15 


0.33 


0.46 


0.74 


0.068 


0.15 


0.36 


0.53 


0.081 


0.27 


0.83 


1.3 


A0.9N50 


5.9 


0.077 


0.17 


0.3 


0.47 


0.59 


0.075 


0.16 


0.33 


0.62 


0.09 


0.3 


0.85 


1.3 


A0.9N100 


5.3 


0.07 


0.16 


0.29 


0.47 


0.51 


0.068 


0.16 


0.31 


0.44 


0.083 


0.29 


0.86 


1.3 


A0.9N200 


5.2 


0.068 


0.15 


0.3 


0.44 


0.49 


0.066 


0.15 


0.39 


0.55 


0.08 


0.29 


0.87 


1.3 


A0.99N100 


8.2 


0.11 


0.19 


0.33 


0.45 


0.59 


0.1 


0.2 


0.38 


0.56 


0.11 


0.34 


0.88 


1.3 



these are unbound outflows with roughly constant rj by the mea- 
sured radius. BH and jet efficiencies are shown decomposed into 
EM, MAKE, PAKE, and EN terms. The radiative efficiency (rj m ) 
for the Novikov-Thorne (NT) model is shown for comparison. 

Many models with \a/M\ > 0.9 show greater than 100% effi- 
ciency (and up to about 300% for the 3D models) for the BH en- 
ergy extraction and jet at larger radii. This is much higher than the 
efficiencies of order the NT efficiency seen in other 3D GRMHD 
simulations that start off with limited poloidal magnetic flux (e.g. 
|Hawley & Krolik|2006> . In the thick disk models, the MAKE ef- 
ficiency is negative even though the particle term (?7^ AKE Po«,) 

is positive. This negative ?/^ AKE is due to a high specific enthalpy 
(i.e. rj™ oc (u g + p g )/po > 1). The negative ^ AKE means the BH 
is accreting thermo-kinetically unbound material as can occur in a 
Bondi flows or does occur in ADAFs (Naray aiT&~Yi|1994| >. For the 
thick disk models, part of the negative 77^ AKE is due the initial torus 
being marginally unbound, which contributes 26% to ?7 H IAKE . 

The dependence upon initial /? is relatively weak for the thick 
disk models, fi varies by factors of 2-3, while the efficiency 77 oc B 2 r 
varies by much less than factors of 2-3 expected if the initial 



fi oc \jB 2 r completely controlled the flux on the BH. T oc B r (dis- 
cussed later) shows even less dispersion vs. initial (/3)~ 1/2 . Similar 
insensitivity to initial fi is seen in Igumenshc hev et aL] ( |2003) . 

Also, notice there are very similar results between the flipping 
and non-flipping poloidal field, which shows that the flipping mod- 
els have plenty of constant polarity flux unlike the MB09D model. 

The 2D axisymmetric simulations (for all thick disk poloidal 
cases, but only showing A0.94BfN40c5* in tables) show relatively 
higher efficiencies than otherwise similar 3D models. Mass is not 
accreted except during infrequent penetrations of the magnetic bar- 
rier via magnetic reconnection. 

Despite the presence of well-ordered poloidal field near the 
BH, we see no evidence for significant energy extraction by any 
ergospheric type disk threaded by magnetic flux that would convert 
spin energy to MA energy in the disk and then to EM energy out in 
the magnetic field (Punsly & Coroniti 1990) unlike suggested to be 
present sometimes in other simulations ( [Punsly et al.]2 009l. Such 
an effect would lead to rj^ AKE » (i.e. outgoing MAKE energy 
on the horizon). For this effect to be dynamically important, one 
should find rj™ AKE > 77™, while in all cases we find 7/Jf AKE « rj™ 
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Table 4. Rest-Mass Accretion and Ejection Rates 



ModelName 


M H 


Mini - 


Mm,o - Mn 


M s 










V A II IIM'M Ilk 


48 


20 


250 


6.9 


6.8 


17 


19 


250 


AU.y4BiJN lUOcl 


37 


16 


170 


6.2 


4.6 


14 


15 


170 


AU.y4BlIN lUUCz 


38 


i 

lj 




A A 

0.0 


5.1 


20 


13 


210 


AU.y4BiJN luUcJ 


39 


17 


200 


6.5 


6.3 


20 


15 


190 


A0.94BfN100c4 


38 


18 


180 


5.4 


6.4 


16 


16 


190 


A0.94BfN40c5* 


4 


15 


280 


1.7 


4.3 


11 


12 


280 


A-0.94BfN40HR 


49 


21 


250 


6.6 


6.5 


18 


20 


240 


A-U.y4BtiN.5U 


v8 


20 


Z50 


A D 
0.8 


o.O 


Zo 


17 


250 


A-U.5BrN30 


1 10 


6.8 


220 


2e-6 


0.53 


7.2 


3.8 


220 


AU.UB1JN 1U 


i An 
lOU 


3.0 


oz 





U.UU4y 


0.033 


5. 1 


100 


AU.jBiN JO 


1 10 


7.6 


140 





0.2 


3.6 


4.3 


140 


AU.y4BiJN JU 


41 


1 8 


160 


7.6 


6.4 


21 


15 


150 


A0.94BfN30r 


41 


17 


170 


5.7 


6.4 


19 


15 


160 


a r\ ct at) xt1 r\r\ 

AU.94BpNlUU 


46 


21 


210 


9.6 


7.2 


24 


19 


200 


A-0.94BtN10 


280 


1.9 


230 





0.0044 


1.7e-6 


1 


230 


A-U.5BtJN 10 


93 


2.2 


84 





0.017 


0.000/ 


2 


86 


AU.UBtJN 1U 


180 


5 







0.089 


o nnnoi 
U.UUUoJ 


A 1 

4. 1 


220 


AU.5BtJN 1U 


200 


5.5 


1 80 





0.29 


0.043 


4.6 


180 


AU.y4BtJN 1U 


190 


16 


420 


0.00095 


0.46 


0.16 


14 


410 


k ft ft A D+VT1 AOD 

Au.y4BtiN 10HK 


330 


-0.83 


410 


0.0085 


1 .6 


0.59 


9.6 


440 


a in nnT~\ 

MBUyD 


0.082 


0.62 


-0.016 


0.0002 


0.001 1 


0.015 


0.77 


0.35 


MBUyQ 


0.94 


0.65 


0.13 


3.1e-5 


0.02 


0.061 


0.74 


1 .9 


A-0.9JN 100 


1 8 


1 .6 


18 


0.32 


0.76 


7.2 


1 .4 


20 


A r\ cm i nn 

A-0.5JN 100 


15 


1.4 


2 1 


0.15 


0.54 


6.8 


1 .3 


25 


A-0.2N100 


15 


1.3 


19 


0.041 


0.31 


5.5 


1.1 


23 


A0.0N100 


17 


1.3 


14 


0.015 


0.4 


6.7 


1.2 


19 


A0.1N100 


14 


1.3 


14 


0.028 


0.36 


5.7 


1.2 


19 


A0.2N100 


15 


1.3 


14 


0.063 


0.35 


5.8 


1.2 


19 


A0.5N100 


13 


1.8 


16 


0.27 


0.72 


g 


1.7 


21 


A0.9N25 


13 


4.3 


20 


0.8 


2.2 


15 


3.7 


22 


A0.9N50 


16 


5.6 


23 


0.88 


3 


16 


4.9 


25 


A0.9N100 


11 


3.9 


21 


0.61 


2 


11 


3.5 


24 


A0.9N200 


7.5 


2.8 


16 


0.5 


1.5 


8.6 


2.4 


17 


A0.99N100 


9.8 


3.7 


14 


0.4 


2.2 


11 


3.3 


16 



except MB09 models give r/^ AKE > r]™. However, in the MB09 
models all of that horizon MAKE energy is dissipated in the disk. 
So in all cases, the jet power is dominated by the EM term, i.e. 
77™ ~ 7/j . In summary, we find that the energy reaching large radii 
is dominated by the EM power produced via the magnetic flux pen- 
etrating the horizon as in the BZ mechanism. 



6.4 Energy Efficiency of Magnetized and Entire Winds 

Table [6] is similar to Table [5] but the magnetized wind ("mw") 
and entire wind ("w") efficiencies (both EM and MAKE decom- 
positions) are shown as evaluated at r, and r given earlier. The 
BH and jet are dominated by EM power, the magnetized wind has 
EM power similar to MAKE power, and the wind is dominated 
by MAKE power (especially at larger radii). The magnetized wind, 
which can reach large radii, is fairly efficient in the poloidal models. 
However, the efficiency is much less than the BH and jet efficiency. 
So the BH dominates the disk's magnetized wind, unlike in some 
models ( |Ghosh & Abramowicz|1997l|Livio et al.|1999> . 



6.5 Angular Momentum Flux for Hole, Jet, and Winds 

A table (not shown) similar to Table [5] except for the specific an- 
gular momentum flux (as computed via the third row of Eq. (14)), 
shows that the poloidal field models have net extraction of angular 
momentum from the BH. The jet carries most of the angular mo- 
mentum and most of that is in EM form. The toroidal field models 
have small angular momentum flux because no steady magnetized 
winds or jets emerge. A table (not shown) similar to Table [6] ex- 
cept for the specific angular momentum flux, shows that the wind's 
angular momentum flux follows our prior discussions for the effi- 
ciency of the winds, except that the EM term tends to dominate the 
MA term for the magnetized wind. 



6.6 Spin-Up Parameter 

Table [7] is similar to Table [5] except the BH spin-up parameter is 



computed via Eq. 1I81, where 



+ s" 



and s K 



s . One can recover the specific angular momentum flux from this 
spin-up parameter and efficiency given in Table [5] The thick disk 
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Table 5. Percent Energy Efficiency: BH, Jet, Winds, and NT 



ModelName 


m 


„™ 


rt MAKE 
''H 


n PAKE 




1) 




n MAKE 
'j 


7?mw,o 


7w,o 


TOT 


A A ti IIX'M lit 


235 


264 


-29.2 


63.4 


-92.5 


242 


228 


13.5 


27.9 


-6.48 


17.9 


AU.y4BtN lUUcl 


323 


379 


-56.1 


83.1 


-139 


333 


322 


1 1 


26.4 


-7.21 


17.9 


A0.y4BiJN 1 UUci 


298 


340 


-42.3 


74.2 


-116 


304 


293 


1 1.3 


34.1 


-4.45 


17.9 


AU.y4.BilN lUOci 


295 


337 


-41.7 


74.3 


-116 


299 


28S 


10.7 


35.5 


-2.95 


17.9 


A().94BfN100c4 


323 


368 


-44.1 


73.2 


-117 


323 


313 


9.44 


31.6 


4.21 


17.9 


A0.94BfN40c5* 


1170 


1240 


-65.5 


86.4 


-152 


945 


925 


19.9 


109 


205 


17.9 


A-0.94BfN40HR 


244 


274 


-30.4 


64.3 


-94.7 


251 


238 


12.1 


27.5 


-7.01 


3.85 


A-U.y4DiJN3U 


87.8 


122 


-34.3 


73.7 


-108 


101 


95.3 


5.38 


15.6 


-12.6 


3.85 


A-U.5BiJN JO 


-13.2 


39.4 


-52.7 


100 


-153 


-1.51 


3.27 


-4.78 


1.15 


-13.7 


4.51 


a n nT> Tin 

AU.UB1JN 1U 


- 17 


-0.1 14 


-16.9 


24.1 


-41 











-5.94 


-12.2 


5.72 


AU.jBiJN jU 


-17.1 


38.8 


-55.9 


85.8 


-142 


-0.796 


5.59 


-6.38 


-6.9 


-18 


8.21 


AU.y4nilN.3U 


353 


420 


-67.1 


88.6 


-156 


361 


348 


12.8 


32.6 


-1 1.4 


17.9 


A().94BfN30r 


323 


386 


-62.7 


86.3 


-149 


334 


326 


8.37 


26.1 


-14.4 


17.9 


AU.y4npJN 1UU 


238 


29 1 


-52.8 


83.8 


-137 


24 1 


220 


20.5 


4U.O 


1 .94 


17.9 


A-0.94BtN10 


-17.8 


-0.703 


-17.1 


24.7 


-41.8 











-9e-7 


-17.6 


3.85 


A-0.5BtN10 


-13.6 


-0.287 


-13.3 


19.8 


-33.1 











-7e-5 


-13.2 


4.51 


AO.OBtN 10 


-18.7 


-0.077 


-18.7 


21.9 


-40.6 











-0.006 


-18.7 


5.72 


A0.5BtN 10 


-19.2 


-0.246 


-19 


26.8 


-45.8 











-0.018 


- 19 


8.21 


AU.y4BtJN 1U 


-22.8 


-0. 162 


-22.7 


33.2 




U.UUU8 


U.UUU4 


U.UUU4 


0.033 


-23.4 


17.9 


AU.y4BtiN 1UHK 


-28.2 


-0.569 


-27.6 


8.71 


-36.3 


0.007 


0.006 


0.002 


0.032 


-19.5 


17.9 


A 1 1 1 1 1 1 1 1 '\ 

MB09D 


14.4 


3.08 


1 1.4 


24.9 


-13.5 


2.96 


2.7 


0.269 


0.485 


-0.739 


16.7 


MB09Q 


5.19 


0.199 


4.99 


22.9 


-17.9 


0.008 


0.004 


0.003 


0.431 


0.576 


17.9 


A-U.yjN 1UU 


32.3 


39.8 


-7.51 


74.6 


-82.2 


23.4 


20.8 


z.ol 


8.49 


9.59 


3.9 


A n C\T 1 (111 

A-0.5JN 100 


17.2 


1 1 .4 


5.77 


41.7 


-36 


9.4 


8.69 


0.706 


7.18 


8.03 


4.51 


A-0.2N100 


6.64 


0.906 


5.73 


34.5 


-28.8 


1.12 


1.04 


0.079 


4.78 


5.38 


5.15 


A0.0N10O 


4.51 


-0.884 


5.4 


33.7 


-28.3 


0.052 


0.035 


0.017 


3.78 


4.23 


5.72 


A0.1N100 


5.04 


-0.242 


5.29 


32.4 


-27.1 


0.342 


0.292 


0.051 


4.16 


4.5 


6.06 


A0.2N100 


7.34 


1.89 


5.46 


34.9 


-29.4 


1.73 


1.64 


0.088 


5.04 


5.38 


6.46 


A0.5N100 


30.8 


24.3 


6.48 


38.1 


-31.6 


19.6 


18.3 


1.25 


10.3 


10.8 


8.21 


A0.9N25 


113 


112 


1.46 


53.3 


-51.8 


88.8 


80.3 


8.55 


24.1 


25.7 


15.6 


A0.9N50 


102 


99.3 


2.91 


51.3 


-48.3 


77.7 


70.3 


7.47 


24.4 


24.8 


15.6 


A0.9N100 


101 


97.2 


4.19 


50 


-45.8 


77.6 


70.1 


7.44 


24.1 


24.7 


15.6 


A0.9N200 


121 


117 


4.28 


52 


-47.7 


96.1 


87.5 


8.67 


25.1 


26 


15.6 


A0.99N100 


130 


127 


2.89 


64.4 


-61.5 


103 


95.5 


7.41 


26.4 


27 


26.4 



poloidal models have a spin-up parameter that is extremely negative 
relative to a/M, unlike the NT thin disk and MB09 models. The 
BH is always rapidly spinning down in absolute spin magnitude 
for all models except the MB09D and the A0.5BtN10 toroidal field 
model. In all cases, the a/M = models are spinning up, although 
quite weakly with the model AO.OBfNIO as compared to models 
AO.OBtNIO and A0.0N100. In summary, the thermal pressure and 
magnetic field are capable of greatly decreasing the BH's spin. 



6.7 a Viscosity and Suppression of the MRI 

Table [8] shows the a viscosity parameters computed via Eq. $29\ 
using various averaging as described below Eq. \29\ in section |3J| 
where a a includes the flow with b 2 /po < 1 and at focuses on the 
highest density parts of the disk flow by using averaging weight 
w = p. a is averaged over r = 12-20r g for all models except 
MB09 models that are averaged from r = 10-12r g . If a compo- 
nent's contribution to a is less than 10%, then we set that compo- 
nent to ~ in the table. Also provided is Q p _ COI (p = n, I, m) com- 



puted via Eq. 1 44 1 for quantities p and b 2 atr = 8r g for the "Disk" 
component. We also show 2a.mri, 2<*,mri> an d Srf.MRi at r dcdcn and 



given in Table 



10 



as computed via Eq. 1 43 i and Eq. 1 33 1. The 

Qnlm .con 2e,MRi! and G^.mri) 



The radius (r s 



i/,MRI = 1 /2 



first occur- 



convergence quality measures (a n 
are discussed in section ( 
rence of 5rf_MRi =1/2 outside 4r g and up to a one inflow radius) 
within which the MRI is suppressed is provided, where "-" indi- 
cates that 5rf,MRi > 1/2 implying no strong MRI suppression. S^mri 
and Srf.weak.MRi (focusing more on heavy disk) give similar results. 

The radially-averaged a depends upon the accumulated mag- 
netic flux, disk thickness, and BH spin. Within the high-density 
disk, the thick disk poloidal field models have \ai,\ ~ 0.02, thick 
disk toroidal field models have a b ~ 0.05-0.07, and thinner disk 
poloidal field models have a b ~ 0.05-0.2. Our thick toroidal field 
models have higher a than found in other works for which the MRI 
is not suppressed, such as a ~ 0.025 found by Beckwit h et al.| 
( |2011} . Also, (not shown) a a > a b , so the less dense and higher 
magnetized corona has larger a as expected. MB09D, that has de- 
caying turbulence, has a b ~ 0.01. 

a and a cff are often quite different for our new thick/thinner 
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Table 6. Percent Energy Efficiency: Magnetized Wind and Entire Wind 



ModelName 




EM 
'mw.i 


„MAKE 
'mw,i 


*7mw,o 


EM 
/mw.o 


/mw,o 


'Av.i 


EM 

'W.I 


„MAKE 
'w,i 




EM 
/w,o 


/w,o 


A A <i i in 

AU.y4BiJN4U 


7.49 


5.9 


1.59 


27.9 


13 


14.9 


-20.1 


6.66 


-26.8 


-6.48 


13.5 


-20 


A0.y4BilN lOOcl 


8.89 


6.5 


2.39 


26.4 


12.3 


14.2 


-19.9 


7.24 


-27.1 


-7.21 


8.16 


-15.4 


A0.y4BilN lOOci 


9.39 


7.36 


2.03 


34.1 


14.6 


19.5 


-18.3 


8.14 


-26.5 


-4.45 


9.99 


-14.4 


A0.y4BlJN lUUcJ 


11.6 


8.07 


3.57 


35.5 


15.4 


20.1 


-17.8 


8.82 


-26.6 


-2.95 


1 1.6 


-14.5 


A0.94BfN100c4 


11.9 


6.98 


4.92 


31.6 


14.8 


16.8 


-20 


7.77 


-27.8 


4.21 


12.3 


-8.08 


A0.94BfN40c5* 


-4.03 


15.4 


-19.5 


109 


43.6 


64.9 


-33.6 


15.9 


-49.5 


205 


42.7 


163 


A-0.94BfN40HR 


8.47 


5.57 


2.9 


27.5 


12.6 


14.9 


-20.4 


6.44 


-26.9 


-7.01 


13.1 


-20.1 


A-U.y4BtJN JU 


3.36 


3.71 


-0.346 


IC ^, 

1 J.O 


7.55 


8.01 


1 O 


A l\A 

4.U4 


-23 


-12.0 


7.62 


-20.2 


A-0.5BiJN JO 


-1.19 


2. 1 


-3.29 


1.15 


6.48 


-5.33 


-13.2 


2.19 


-15.4 


-13.7 


6.06 


-19.7 


A0.0B1JN10 


-13.3 


0.068 


-13.4 


-5.94 


0.087 


-6.03 


-16.8 


0.087 


-16.9 


-12.2 


0.26 


-12.4 


AO.jBlJN JO 


-3.82 


1.54 


-5.3 / 


-o.y 


3.92 


-10.8 


-14.8 


1 .42 


-io.2 


-18 


A nn 


T"> 1 
-11. 1 


A l~\ n AT> 4Tv T "5 

A0.y4BlJN JO 


11.2 


7.26 


3.9 


32.6 


15.4 


17.2 


-19.7 


7.89 


-27.6 


-1 1.4 


16.7 


-28.1 


A0.94BfN30r 


9.03 


6.62 


2.41 


26.1 


11.9 


14.2 


-21 


7.14 


-28.1 


-14.4 


13.3 


-21.1 


AU.y4BpJN 1UU 


10.8 


7.24 


3.52 


40.6 


18.4 


22.2 


1 ^ n 
-16.9 


8.23 


-25.1 


1 .94 


20.6 


-18.7 


A-0.94BtN10 


-0.164 


-0.017 


-0.147 


-9e-7 


-2e-7 


-7e-7 


-17.6 


-0.148 


-17.4 


-17.6 


0.019 


-17.7 


A-0.5BtJN 10 


-0.392 


-0.014 


-0.378 


-7e-5 


-3e-5 


-4e-5 


-13.5 


-0.081 


-13.4 


-13.2 


0.04 


-13.2 


AO.OBtJN 10 


-0.457 


0.026 


-0.483 


-0.006 


-9e-5 


-0.006 


-18.7 


0.054 


-18.8 


-18.7 


0.057 


-18.7 


AU.jBtJN 10 


-1.36 


0.077 


-1.44 


-0.018 


0.003 


-0.022 


-ly 


0.098 


i n 1 
-iy.i 


-19 


0.072 


1 n 1 

-iy.i 


A0.y4BtN10 


-0.413 


0.09 


-0.503 


0.033 


0.012 


0.021 


-22.9 


0.109 


-23 


-23.4 


0.113 


-23.5 


AU.V4rStiN luHK 


-1.66 


0.162 


-1.83 


0.032 


0.04 


-0.008 


-23. 1 


0.187 


-23.3 


1 n ^ 
-19.5 


0.177 


-19.7 


MBuyiJ 


3. ID 


3.42 


-0.271 


0.485 


1 .28 


-0.792 


5.99 


o.lo 


-0.184 


-0.739 


1.30 


-2.09 


MBuyy 


2.49 


1.65 


0.839 


0.431 


0.985 


-0.555 


4.84 


2.17 


2.67 


0.576 


1.14 


-0.567 


A-0.yN100 


3.97 


2.26 


1.71 


8.49 


3.72 


4.77 


5.41 


2.28 


3.12 


9.59 


3.85 


5.75 


A r\ CM 1 Art 

A-0.5N 100 


3.57 


2.05 


1.52 


7.18 


3.39 


3.78 


5.3 


2.07 


3.23 


8.03 


3.51 


4.52 


A-0.2N100 


2.85 


1.54 


1.32 


4.78 


2.52 


2.26 


4.82 


1.56 


3.26 


5.38 


2.61 


2.77 


A0.0N100 


2.88 


1.67 


1.21 


3.78 


2.37 


1.41 


4.5 


1.7 


2.79 


4.23 


2.45 


1.78 


A0.1N100 


2.97 


1.73 


1.24 


4.16 


2.78 


1.38 


4.56 


1.75 


2.81 


4.5 


2.82 


1.69 


A0.2N100 


3.42 


2.09 


1.32 


5.04 


3.42 


1.63 


4.94 


2.1 


2.83 


5.38 


3.49 


1.89 


A0.5N100 


5.77 


4.3 


1.47 


10.3 


6.85 


3.46 


7.45 


4.33 


3.12 


10.8 


6.95 


3.83 


A0.9N25 


9.8 


7.01 


2.79 


24.1 


12.2 


11.9 


11.7 


7.09 


4.66 


25.7 


12.3 


13.4 


A0.9N50 


10.6 


7.86 


2.78 


24.4 


13.1 


11.3 


13 


7.94 


5.03 


24.8 


13.1 


11.7 


A0.9N100 


11 


8.21 


2.84 


24.1 


13.3 


10.8 


13.4 


8.27 


5.13 


24.7 


13.4 


11.3 


A0.9N200 


11.4 


9.09 


2.26 


25.1 


14.3 


10.8 


13 


9.22 


3.73 


26 


14.5 


11.5 


A0.99N100 


11.8 


8.69 


3.08 


26.4 


13.8 


12.5 


13.6 


8.79 


4.85 


27 


14 


13 



poloidal field models and thick toroidal field models, which shows 
that local viscous a-disk theory is inaccurate. For the poloidal field 
models, o>disk theory fails for r < 10r g due to disk compression 
(leading to smaller ff 1 than a hydrostatic solution would have) by 
magnetic flux. Even at larger radii within the inflow equilibrium 
region where ff 1 is not compressed much by magnetic flux, cf-disk 
theory still does poorly. Using a b oc p b leads to somewhat better 
agreement for some of the poloidal field models. The overall lack 
of agreement shows that large-scale magnetic torques through mag- 
netic confinement (rather than direct magnetic torques within the 
heavy inflow) play an important role. Convection may also play an 
important role. For the toroidal field models, using at oc p tot gives a 
radial scaling for v vix similar to v r but too small by a factor of ~ 10, 
while a b oc p b gives far too large v vix compared to v r . Summariz- 
ing, in all our new models, using p tot leads to too small v r = u v ; sc 
despite often having a reasonable radial scaling. Our poloidal and 
toroidal field simulations have a eB ~ 0.2-1 (which indicates the 
actual effective viscous timescale) that is large enough to be con- 
sistent with observations (King et al. 2007 ). 

By contrast, the same table diagnostics for the multi-field loop 



thin (ff 1 ~ 0.07) disk model A0HR07 in |Penna et aT] ( |2010| > give 
a b ~ a b.cff ~ 0.04 for r = 7-9r g , so a-disk theory works quite well 
for those multi-loop field thin disk models. Also, for ff' ~ 0.07, 
|Beckwith et al.H201 if found a ~ 0.025 and a eff ~ 0.1 (from their 
fig. 8 around r ~ I0r g where G ~ 2), so a-disk theory is holding (at 
best) marginally well. The tables show that our older MB09 models 
show reasonable agreement with o/-disk theory. 

Ob,M2 dominates the local stress contribution to at in many 
cases. Note that averaging |a- mag | gives values of ~ 0.4-0.6 for all 
our simulations (including 2D models), but that neglects the direc- 
tion of angular momentum transport. The PA (the Reynolds stress) 
term can be computed as qtpa ~ a - aM2 from the table since other 
terms are small. Interestingly, at, a bi Mi, and «A,mag are negative for 
a/M < for the thick disk poloidal field models, and such models 
even have v Iot the same sign as a/M (i.e. flow reversal) in the heavy 
disk inflow even out to r ~ 40r g (e.g. in model A-0.94BfN30). This 
flow reversal behavior was validated by restarting the A0.94BfN40 
model at t = 8000r s /c with a/M = -0.9375 that produced model 
A-0.94BfN40HR. (Compared to A-0.94BfN30, A-0.94BfN40HR 
shows much less variance in time-dependent quantities - proba- 



© 2012 RAS, MNRAS 000,[T]{35] 



Magnetically Choked Accretion Flows 25 



Table 7. Spin-Up Parameter: BH, Jet, Winds, and NT 



ModelName 




..EM 
' S H 


( ,MA 
J H 


pEA 
' S H 


..EN 
4 H 


s ) 


v 


j 


Wo 




-VNT 


it A li IIX'M lit 


-23.5 


-21.1 


-2.37 


-0.609 


-1.76 


-19.5 


-16.8 


-2.65 


-3. 1 


-4.03 


0.41 1 


AU.y4BiJN lUUcl 


-28 


-25 


-3.01 


-0.27 


-2.74 


-22.6 


-20.3 


-2.31 


-3.47 


-5.67 


0.41 1 


AU.y4BiJN lUUci 


-26.6 


-23.6 


-2.95 


-0.478 


-2.47 


-21.1 


-18.6 


-2.44 


-3.45 


-2.92 


0.41 1 


A A A /I ~F» -TNT 1 aa „ 1 

AU.y4rSlJN 1U0c3 


-25.7 


-22.7 


-2.97 


-0.484 


-2.48 


-20.2 


-17.8 


-2.43 


-4.23 


-4.44 


0.41 1 


A0.94BfN100c4 


-26.8 


-24.1 


-2.75 


-0.468 


-2.29 


-20.8 


-18.4 


-2.34 


-3.11 


-6 


0.411 


A0.94BfN40c5* 


-111 


-105 


-6.41 


-0.747 


-5.66 


-107 


-103 


-4.75 


-11.6 


-67.1 


0.411 


A-0.94BfN40HR 


23.5 


21.1 


2.31 


0.573 


1.74 


19.4 


16.9 


2.49 


3.19 


7.76 


6 


A-U.y4BiJN30 


11.4 


9.54 


1.89 


0.264 


1.63 


8.81 


6.67 


2.14 


2.55 


5.01 


6 


a a c n fMin 

A-U.5B1N 30 


12.7 


1 1.5 


1.23 


-0.144 


1.37 


5.62 


4.91 


0.714 


6.97 


8.5 


4.84 


A A AFJ -CM 1 A 

AO.OBilN 1U 


0.013 


0.032 


-0.019 


-0.009 


-0.009 


-0 


-0 


-0 


0.061 


0.404 


3.46 


AU.jbtJN 3U 


-11.5 


-10.2 


-1.3 


-0.005 


-1.3 


-5.15 


-4.55 


-0.599 


-D.j4 


-8.05 


1 .98 


A A A /I T» -CNTf A 

A0.y4rSlJN30 


-30.6 


-27.5 


-3.08 


-0.162 


-2.92 


-25 


-22.3 


-2.62 


-4.06 


-10.3 


0.41 1 


A0.94BfN30r 


-28.9 


-25.9 


-3 


-0.205 


-2.79 


-23.1 


-20.7 


-2.45 


-3.63 


-8.54 


0.411 


Afl (lylD «\i 1 a a 
A(J.y4i3pJN 1UU 


-25 


-22.3 


-2.72 


-0.214 


-2.5 


-21.3 


-18.7 


-2.56 


-5.07 


-5.97 


0.41 1 


A-0.94BtN10 


3.29 


0.205 


3.09 


1.98 


1.11 


1.88 





1.88 


1.88 


3.2 


6 


A II CI') *M 1 fl 

A-0.5BtN10 


2.39 


0.103 


2.29 


1.61 


0.678 


1 





1 


1 


2.18 


4.84 


AU.UBtNlU 


1.15 


0.004 


1.15 


0.751 


0.396 


-0 


-0 


-0 


0.0004 


1.06 


3.46 


A n en *M 1 1~\ 

A0.5BtN10 


-0.639 


-0.297 


-0.342 


-0.208 


-0.134 


-1 


-0 


-1 


-1 


-0.793 


1.98 


A0.y4BtJN 10 


-1.53 


-0.232 


-1.3 


-0.716 


-0.584 


-1.88 


n nn/vi 

-0.0003 


-1.87 


-1.88 


-1.63 


0.41 1 


A A Ik 1 l>i\[ 1 HI 1 11 

AU.y4BtJN 1UHK 


-2.04 


-0.316 


-1.72 


-1.38 


-0.336 


-1.88 


-0.002 


-1.87 


-1.88 


-2.33 


0.41 1 


A in nun 


-0.204 


-0.069 


-0.135 


-0.039 


-0.096 


-1.99 


-0.137 


-1.85 


-8.63 


-15.3 


0.492 


md nnn 
MBUytj 


0.103 


-0.046 


0. 149 


0.132 


0.017 


-1.88 


-0.0005 


-1.87 


-5.39 


-7.63 


0.41 1 


a n nM i nil 
A-U.yjN 1UU 


5.92 


3.17 


2.75 


0.896 


1,00 


3.03 


1.15 


1.88 


-0. 15 1 


3.52 


C II 


A 11 CM I nn 

A-0.5JN 100 


5.28 


3.63 


1.65 


1.03 


0.625 


2.61 


1.54 


1.06 


-0.819 


1.52 


4.84 


A-0.2N100 


3.4 


2.19 


1.21 


0.816 


0.399 


1.17 


0.765 


0.404 


-0.761 


0.879 


4.02 


A0.0N100 


1.03 


0.004 


1.03 


0.726 


0.304 


-0.045 


-0.045 


-0.0004 


-3.19 


-0.6 


3.46 


A0.1N100 


-0.546 


-1.03 


0.481 


0.339 


0.141 


-0.649 


-0.448 


-0.201 


-3.91 


-2.35 


3.18 


A0.2N100 


-1.53 


-1.98 


0.451 


0.323 


0.128 


-1.37 


-0.954 


-0.418 


-4.28 


-2.41 


2.89 


A0.5N100 


-5.05 


-5.15 


0.099 


0.096 


0.002 


-4.43 


-3.27 


-1.16 


-6.39 


-4 


1.98 


A0.yN25 


-8.5 


-7.31 


-1.18 


-0.511 


-0.674 


-6.5 


-4.33 


-2.17 


-7.77 


-4.49 


0.58 


A0.9N50 


-7.46 


-6.44 


-1.01 


-0.45 


-0.565 


-6.05 


-3.9 


-2.15 


-8.34 


-4.72 


0.58 


A0.9N100 


-7.28 


-6.35 


-0.925 


-0.429 


-0.496 


-6.05 


-3.92 


-2.13 


-8.74 


-5.41 


0.58 


A0.9N200 


-8.64 


-7.57 


-1.07 


-0.477 


-0.591 


-7.24 


-4.98 


-2.26 


-10.1 


-4.42 


0.58 


A0.99N100 


-6.8 


-5.57 


-1.22 


-0.394 


-0.828 


-5.51 


-3.34 


-2.16 


-9.7 


-4.57 


0.111 



bly due to lower-resolution models allowing some of the oppo- 
site polarity magnetic flux to reach the BH causing, e.g., the time- 
averaged ?7h to be smaller.) Recall that o- mag oc —b r b$, which is 
just a term related to the Poynting flux. Evidently, BH angular 
momentum is being dumped even into the heavy disk. Yet, in all 
cases, mass continues to accrete due to large-scale stresses by mag- 
netic flux and also possibly due to convection. Interestingly, the 
thick disk toroidal field models are dominated by positive a b px 
(Reynolds stress). The thick disk a/M = poloidal field model 
AO.OBfNIO has non-negligible negative a bPA . Also, some of the 
thinner disk poloidal field models have (negative) Reynolds stress 
that is up to 50% of lo^jvol- ff&.EN contributes little to art, except for 
the thick disk toroidal field models for which a b w < a bM2 < both at 
about 20% contribution to the total a. 

The MRI is not suppressed (Sj.mri ^ 1) for thick disks at 
low spin, flows with initially toroidally-dominated field, or MB09 
type limited poloidal flux or quadrupolar field models. However, 
for poloidal field models where magnetic flux saturates, the MRI 
is suppressed within t's dMRI =i/2 ~ 20-200r g for many thick disk 
models and within r SjMSI=i / 2 ~ 20r g for the thinner disk models 



(the radius depends upon the initial /? and duration of simulation). 
The MRI is suppressed in the AO.OBfNIO model as due to the small 
initial /? and much of the magnetic field staying in the dense flow 
(instead of a jet) where S^mri is measured. 

The suppression of the MRI is due to both the field strength 
and the angular velocity. The thinner disk models have accumu- 
lated much magnetic flux near the BH as most evident from Fig- 
ure [T9]( similar to Figure |6j showing the time-averaged flow-field 
for A0.99N100. This suggests that thinner disks accrete (or hold 
onto) magnetic flux as easily as thicker disks. However, the thick 
flows are more sub-Keplerian, as evident when computing S^mri 
using Q rot — > £1k- Then, the thinner models would still have 
rs dMRI =i/2 ~ 20r s because they are somewhat Keplerian, while 
the thick models would have ^s rfMRI =i/2 ~ 8r g . The MB09Q and 
toroidal field models have 5rf_MRi ~ few because poloidal flux build- 
up is competing against its destruction since there is no net flux 
except in small patches. Our high resolution toroidal field model 
shows smaller 5 f ;,MRi, so even higher resolutions might lead to 
SjMRi ~ 0-3 as in the poloidal field models. MB09D probably has 
Srf,MRi ~ 8 simply because turbulence is under-resolved. 
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Table 8. Viscosities, Grid Cells per Correlation lengths and MRI Wavelengths, MRI Wavelengths per full Disk Height, and Radii for MRI Suppression 



ModelName 


»6,eff 




«Z),M2 




&.,cor, 


2l,cor, 


&>,cor, 


ye,MRI.(i,ol 


y0,MRI,{i,o) 


r 

O d,MRI,{i,o) 


1 ii ■ u.weiiK 1 












IPO'* 2 ! 


IPO.* 2 ! 










MRI=l/2 


A A n 1 1J+'Xi in 
AU.y4mJM4U 


0.39 


0.022 


0.02 


0.26 


22, 21 


O, 12 


19, 18 


110, 120 


ion C/in 
joO, j4U 


0.33, 0.34 


> 55, > 55 


Au.y4BtJN lUOcl 


0.35 


0.029 


0.029 


0.36 


13, 13 


4, 7 


12, 1 1 


56, 53 


i i n 

210, 220 


0.29, 0.38 


64, 44 


A0.y4BlJN 1UUCZ 


0.3 


0.029 


0.028 


0.3 


13, 13 


4, 7 


o, o 


A O A 1 

48, ol 


98, 100 


n i q n ia 
0.38, 0.34 


> 44, > 44 


AU.y4BiJN lUOcJ 


0.32 


0.03 


0.026 


0.31 


14, 13 


4, 7 


3, 3 


54, 62 


47, 59 


0.31, 0.33 


47, 42 


A0.94BfN100c4 


0.29 


0.022 


0.018 


0.24 


13, 13 


3, 8 


2,2 


36, 48 


24,31 


0.47, 0.43 


14, 72 


A0.94BfN40c5* 


0.089 


0.0053 


0.0032 


0.051 


22, 19 


3, 22 


1, 1 


1300, 660 


1, 1 


0.025, 0.058 


55,60 


A-0.94BfN40HR 


0.088 


-0.017 


-0.017 


-0.24 


22, 21 


6, 12 


20, 18 


110, 140 


380, 560 


0.36, 0.29 


> 54, > 54 


A-U.y4BlJN JU 


0.12 


-0.04 


-0.031 


-0.37 


15, 13 


-1 "7 

4, / 


10, 9 


A1 AG 

42, 45 


1 on o^n 
loO, 2o0 


C\ A A A AC 

0.44, 0.40 


T 1 A 1 AA 

210, iyo 


A-0.5B1N30 


0.2 


-0.1 1 


-0.11 


-0.65 


13, 13 


5, 6 


15, 16 


35, 23 


inn i c\r\ 

200, 290 


0.5, 0.91 




a n no pnt i n 


0.77 


-0.012 


~ u 


-0.01 1 


15, 15 


12, 8 


13. 13 


110, 120 


110, 170 


n T7 A 1 n 
0.22, o.iy 


. OO v OO 
> OO, > OO 


AO.jBIJNjO 


0.39 


0.12 


0.12 


O.OJ 


14, 13 


3, 7 


16, 19 


47, 18 


240, 2/0 


0.39, 1.2 


18, 18 


a n flvi'Dfvran 
AU.y4BtiNiU 


0.36 


0.024 


0.022 


0.25 


13, 13 


3, 7 


12, 12 


35, 28 


1 1 n i~ir\ 
210, j /O 


0.4b, 0. /j 


13, 23 


A0.94BfN30r 


0.28 


0.024 


0.023 


0.19 


13, 13 


3,6 


12, 12 


31,34 


190, 330 


0.52, 0.6 


11, 13 


a n n /i "d «\t i nn 
AU.y4BpJN 1UU 


0.45 


0.028 


0.024 


0.3 


13, 13 


4. 8 


12,11 


51, 100 


1 on T7fl 
160, 2 /U 


n n A 1 n 
0.J2, o.iy 


1 i a 1 nn 
1J0, 100 


A-0.94BtN10 


0.14 


0.053 


~ 


0.4 


15, 13 


17, 10 


20,6 


19, 10 


40, 32 


1.4,2.7 


-■- 


A n CD4KT1 n 
A-U.jBtJN 1U 


0.17 


0.096 


0.01 


0.38 


15, 13 


18, 10 


20, 7 


23, 12 


48, 40 


1 i i i 
1.2, 2. 1 




AO.OBtJN 10 


0.24 


0.06 


0.016 


0.47 


15, 13 


17, 9 


20, 8 


19, 15 


61, 64 


1.5, 1.7 




A0.5BUN 10 


0.48 


0.067 


0.022 


0.44 


15, 15 


15, 12 


15, 9 


31, 23 


87, 83 


0.9, 1.1 




A0.y4BtJN 10 


0.4 


0.062 


0.00/2 


0.34 


14, 14 


16, 1 1 


19, 10 


19, 16 


58, 63 


1.5, 1.7 




A A III D*KT1 ADD 

Au.94.BtJN 1UHK 


0.75 


0.066 


0.02 


0.38 


28, 25 


29, 24 


25, 18 


71, 50 


170, 1 /0 


0.8, 1 




MBU9JL) 


0.014 


0.01 


0.00/9 


0.23 


18, 14 


4, 4 


4, 3 


3, 3 


4, 4 


7.5, 8.1 




jvibuytj 


0.099 


0.076 


0.054 


0.35 


22, 19 


12, 9 


3, 3 


11,11 


9, 9 


1 1 Ci 

3, 2.y 




a n o\i i r\n 
A-U.yJN 10U 


0.25 


0.034 


0.18 


o.4y 


23, 22 


9, 1 1 


5, 5 


1 10, 35 


22, 21 


a on A O A 
0.2y, O.o4 


22, 20 


\ r\ cm i aa 

A-0.5N 100 


0.51 


0.067 


0.19 


0.48 


32, 31 


9, 9 


4, 3 


170, 33 


25, 22 


n i o r\ o a 

0.18, 0.84 


25, 22 


A-0.2N100 


0.68 


0.12 


0.17 


0.46 


32, 29 


8,7 


3, 3 


140, 22 


23, 21 


0.21, 1.4 


20, 18 


A0.0N100 


0.62 


0.21 


0.18 


0.52 


31, 29 


9, 10 


4, 3 


110,32 


23, 21 


0.28, 1.1 


21, 18 


A0.1N100 


0.67 


0.2 


0.21 


0.49 


31, 31 


9, 10 


3,3 


140, 35 


28, 23 


0.23, 1 


21, 19 


A0.2N100 


0.92 


0.35 


0.23 


0.57 


32, 31 


8, 8 


3, 3 


150, 34 


29, 22 


0.21, 0.94 


24, 23 


A0.5N100 


0.74 


0.24 


0.18 


0.53 


31,28 


3,7 


3, 3 


90, 30 


25,20 


0.32, 1.1 


16, 15 


A0.9N25 


1.2 


0.13 


0.26 


0.59 


29,27 


8, 10 


4,4 


160, 280 


48,45 


0.21,0.16 


> 120, > 120 


A0.9N50 


0.96 


0.16 


0.21 


0.55 


28, 27 


8, 10 


3,4 


130, 88 


30, 23 


0.25, 0.4 


34, 32 


A0.9N100 


0.92 


0.17 


0.21 


0.54 


29, 28 


8, 10 


4,4 


120, 42 


30, 22 


0.27, 0.79 


23,21 


A0.9N200 


0.53 


0.18 


0.15 


0.46 


27,26 


8,9 


3,4 


75,24 


26,24 


0.42, 1.5 


14, 12 


A0.99N100 


0.83 


0.2 


0.2 


0.53 


22,20 


6.8 


7.7 


110, 67 


66, 54 


0.33, 0.56 


26, 20 



6.8 Magnetic Flux 

Table [9] shows the value of Y (computed via the Eq. l |26| >) on the 
horizon (Th), in the inflow-only (u r < 0) regions at r, and r a de- 
fined already (T in ; and Y ino respectively), in the jet (Yj), in the 
magnetized wind at r, and r„ (Y mw ; and Y mw o respectively), and in 
the entire wind at r,- and r„ (Y mw> i and Y mWi0 respectively). Only Y^ 
uses the non-local normalization of M H - 

In the poloidal field models, Y is dominated by the magnetic 
flux threading the BH polar regions as shown by comparing Y H and 
Tj. Also evident is that Y decreases from r„ to r, (shown by Y in o 
and Yj„_j, respectively), but then rises near the horizon. The flux at 
large radii acts a reservoir. At smaller radii some flux gets accreted, 
and on the horizon much of the mass gets drained. This leaves the 
ordered part of the magnetic flux giving high Yh- 

The magnetized wind and entire wind values of Y are normal- 
ized by the local mass accretion rate for the inflow component at 
those radii (i.e. M iu i and M ino ). This shows that the local mass- 
loading can be large in the wind and small in the magnetized wind. 



The entire wind is just part of the overall circularizing flow and 
eventually feeds the pure inflow and feeds the inflow value of Y. 

Also, notice that Y is quite similar between the flipping and 
non-flipping poloidal field, which shows that the flipping models 
have plenty of constant polarity flux unlike the MB09D model. 
This also shows that even with different initial conditions, the flux 
threading the BH saturates in some type of force balance between 
(e.g.) the inner magnetospheric pressure against the exterior disk 
gas+ram pressure. This unstable balance results in time-dependent 
magnetic Rayleigh-Taylor events that regulate the magnetic flux on 
the horizon by carrying magnetic flux back into the disk and wind. 

Table [9] also shows the magnetic flux (computed as in sec- 
tion [3T6J on the horizon normalized in various ways. The values 
Th/T^/ = 0), Th/Tj^ = 0), and Tn/^a = 0) are for normal- 
izations by the extrema (label of 1 for the extrema at smallest radii, 
and so forth) in the magnetic flux in the initial disk. For the flipping 
poloidal field geometry, we show only the first 3 extrema since only 
those are relevant over the entire evolution. A value of zero indi- 
cates no relevant value, as for the non-flipping poloidal field ge- 
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Table 9. Absolute Magnetic Flux per unit: Rest-Mass Fluxes, Initial Magnetic Fluxes, Available Magnetic Fluxes, and BH Magnetic Flux 
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ometries. For the toroidal field models, this just shows the growth 
of the initially small magnetic flux threading the equator. 

The value < F H / < F fl (as other quantities, this is [YhO) /[¥„],], 
in the table) is [the magnetic flux on the horizon] per unit [flux 
available on the BH plus beyond the BH of the same polarity of 
magnetic flux]. This ratio shows how much more flux is available 
to the BH, and a value of ^Vn/^a ~ 1 would mean the magnetic 
flux beyond the BH is only of opposite polarity - so that the BH 
has as much flux as it can get of the same field polarity. The value 
^a/^s is similar to *P H /*P a , except it shows how much magnetic 
flux is available to the BH of any polarity within the stagnation ra- 
dius. The stagnation radius (r,) is defined by where = f ptf/dA&p 
over the full flow within b 2 /p < 1 (so includes the disk inflow and 
entire wind outflows), because the whole flow has not reached in- 
flow equilibrium there. This is roughly where u r = as weighted 
for the whole massive flow. If the density-weighted radial lab-frame 
3-velocity v T = inside this radius, then that radius is used instead 
for r s . Note that y V a and *F, have been time-averaged before form- 
ing a ratio in the table. Lastly, the value IOh/^ihI shows the value 



of the absolute magnetic flux per unit absolute of the extremum 
of the signed magnetic flux, where the ratio itself is time-averaged 
since the measurements are exactly co-spatial. As discussed in sec- 
tion |3.6| this roughly measures the vector spherical harmonic / 
mode (e.g. I = 1 is dipolar). 

For the polarity-flipping poloidal field models, the values of 
^H/^i s£ 1 if the magnetic flux on the BH has already exceeded the 
i-th loop's initial magnetic flux. So, in our polarity flip simulations, 
the table shows that magnetic flux has accumulated and been de- 
stroyed already twice in the simulation with the flipping model set- 
tling on the third extremum (i.e. 3rd field loop). Longer times lead 
to the next polarity being accreted. This shows that ordered flux is 
easily transported through the accretion disk, and large-scale field 
that penetrates the corona (see, e.g., |Rothstein & Lovela ce 2008 1 is 
not required for flux transport. Movies of the field lines (contours of 
A, or streamlines of B') do show that magnetic field lines at higher 
latitudes initially transport more readily than those at the equator 
as seen in more idealized models (Beckwith et al. 20091. However, 
eventually the magnetic flux threading the equator does accrete as 
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-10 



Figure 19. Flow-field (as in Figure[6]for our fiducial thick disk model) for 
model A0.99N100 with a/M = 0.99 and a thinner disk. The thinner disks 
show less magnetic flux threading the BH. The dense inflow accretes onto 
the BH through the ordered magnetic flux via efficient non-axisymmetric 
magnetic Rayleigh-Taylor instabilities. 



well, and in steady state we cannot identify anything special about 
the corona vs. the disk in accretion of ordered magnetic flux. |Beck^| 
|with et a!7| <[2009 ) suggest their simulations may have reached mag- 
netic flux saturation near the BH (i.e. saturated Th and Yh), but 
they also found that most of the magnetic flux moves out to large 
radii rather than accreting. We suggest that this is due to their small 
initial torus leading to significant outflow of magnetic flux (as in 
their fig. 2, fig. 3, and animation). 

The values of ^u/^a and y ¥n/ y ¥ s show that there is plenty 
of magnetic flux available (including of the same polarity) for the 
BH in most models. These are time-averages, while in general each 
polarity loop that is accreted starts off at quite low values of V P H / V P (1 
until nearing the next polarity when v P H / v P n ~ 1 and then the field 
polarity inversion occurs on the horizon. 

The MB09D model, typical of most torus-based MHD simu- 
lations in the literature, was time-averaged over an early period of 



0.5 



-0.5 




2xl0 4 



4xl0 4 6xl0 4 8xl0 4 
t[r g /c] 



Figure 20. Value of'T t H/'I ) H corresponding to the approximate value of I /I 
for the vector spherical harmonic expansion /-mode for model A0.94BtN10 
(similar results for A0.94BtN10HR). This shows that well-ordered dipolar 
(IftH /'I'hI ~ 1, where I'PtH/'l'Hl ~ 1/2 corresponds to quadrupolar) field 
on the BH appears for somewhat sustained durations. In summary, only 
weak relativistic jets are produced in our simulations that start with mostly 
toroidal field, but much higher resolutions might allow this emergent large- 
scale dipolar field to launch persistent lightly-loaded relativistic jets. 



For model A0.94BtN10 (A0.94BtN10HR is similar), Figure [20] 
shows ^h/Oh (approximate 1// for the vector spherical harmonic 
I mode). The large-scale dipolar field is evident in Figure[2T| where 
the large-scale dipolar flux extends out to r ~ 50r g . 

By contrast, Figure [22] shows the typical flow at times when 
I'Pih/ShI * !• The magnetic field is strongly mass-loaded and not 
particularly aligned with any rotational axis. The high-density ma- 
terial inflows from quite random directions. 

During the episodes of large-scale dipolar flux formation in 
the toroidal field models, only weakly powered highly relativistic 
jets emerge due to too low T ~ 3. Higher resolutions, higher T, 
or somewhat thinner disks might promote more persistent jet for- 
mation, which is required by some reconnection-based jet models 



accretion. However, already by t ~ 2000r„/c the >F H /«F -» 1 and JUzdensky & McKinney|20 11 [[McKinney & Uzdensky|2012). 



^h/^Vs — » 1, indicating there is no more available magnetic flux 
with the same polarity (either at all or within the part of the ingoing 
flow, respectively). Only the accretion of another ~ 5x the same 
polarity magnetic flux would have led to a MCAF state and much 
higher rj. Also, some prior MHD simulations may have been run for 
too short of a duration. Depending upon the initial conditions, the 
magnetic flux can accumulate over longer times than other quanti- 
ties, so looking at only energy and angular momentum fluxes can 
be misleading. In short, the choice of initial conditions and short 
duration can limit Th and so the efficiencies. 

The value of lOn/ftHl shows the average approximate vector 
spherical harmonic |/|-mode. All poloidal field models have unity 
as expected, while the MB09Q large-scale quadrupolar field model 
has I ~ 3 as expected given the large-scale field is quadrupolar 
(/ = 2) and there is an equatorial MHD turbulent disk. 

Interestingly, the toroidal models show episodes with well- 
ordered and large-scale dipolar |/| « 1 field on the horizon. 



6.9 Power-Law Fits for Radial Dependence 

The power-law scaling of quantities in the flow is important for 
testing accretion flow models. RIAFs such as ADAFs have p oc 



r~ 3 ' 2 and p g oc r~ 5 ' 2 , while CDAFs have p oc r _1/i , and ADIOSs 
vary over this range. ADAFs have cJv K ~ 0.5-0.7 and cjv$ > 1. 

Viscous simulations and RIAF models show similar radial 
scalings for density, velocity, an d mass accretion rate (Stone et al.j 
|1999| [i gumenshc hev et al.[2 000 , Igumensh chev & Narayan|2002[ 



|McKinney & Gamrnie 



2002 1. For example, flo ws with ff 1 ~ 0.2 



show an inflow M oc r (Hawley & Balbus 20021. For viscous sim- 
ulations, Igumenshchev & Abramowicz||2000| > found p oc r~ [/1 to 
r~ 0J for r = 5/3-4/3, respectively. For relatively thin disk injec- 
tion 3D MHD simulations giving c s /v Iot ~ 0.4, Igumenshchev et al. 



1 2003 i found p oc r 1 and v I0t oc r 1/2 . In 2D MHD simulations with 
Cj/fmt ~ 0-6, Stone & Pringle (2 001[ > (i.e. Run F) found p oc con- 
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Figure 21. Like upper-left panel in Figure [4] and Figure |13| but for ini- 
tially toroidally-dominated magnetic field model A0.94BtN10. This shows 
"spontaneous" formation of large-scale and dipolar magnetic flux near the 
BH at t a 62828r g /c - as also occurs at other times when |T t H /^H I ~ 1 (see 
Figure pO) . The magnetic field lines are shown as thick (thin) black lines 
when they are lightly (heavily) mass-loaded. Only lightly mass-loaded field 
lines would permit ultrarelativistic jets. The magnetic field and low-density 
(blue in plot) wind/jet regions extend out to (not shown) r ~ 50r g . 




Figure 22. Like Figure pT| but typical snapshot at t ss 40164r g /c when 
I'PfH/^Hl *K 1 while the field is not dipolar or large-scale. Magnetic field 
lines are strongly mass-loaded (i.e. field lines shown as thin lines). Field 
lines appear/disappear because (in this slice) they bend in/out of the plotted 
plane. In such thick flows, mass inflow occurs from quite random directions 
with no strong preference for the plane of the disk/BH rotation. 



stant, 



Pg K r 



and v r oc r 1 . However, these MHD simulations 



involved relatively thin disks with 0' p < 0.6. For MHD simulations 
of truly thick (0 ~ 1 or c s /v m » 1) flows, |Pen et al-H2003| found 
p oc r~ 0J2 , p g oc r~ l 5 , and v Iot ~ O.I^k - a sub-Keplerian flow. 

We now consider radial power-law fits (to the form / = 
fo(r/r )") for various quantities. Table [T0| shows some of our re- 
sults. The 5-cr errors in the least squares fits for the power-law index 
are shown as a subscript, where systematic errors (not statistical er- 
rors) dominate the lack of a fit. We also by-eye looked at all fits and 



confirmed that they are reasonable. This is how we tuned the time, 
number of inflow times, and number of <x that best represent some- 
thing about the systematic uncertainties. A "-" is shown if the 5-cr 
error permits n to pass through zero, except for \n\ < 0.25 in which 
case a "-" is shown if the 5-cr error is larger than 0.25. MB09D and 
2D models are not shown because of their small inflow equilibrium 
radii. The resulting power-law can change at larger radii because 
the flow is not strongly self-similar where fits are obtained. 

Table|10|shows fits for "disk quantities": po, p g , and \b\ (where 
l^tyl ~ \b\)- We also consider (not shown) \v r \, \v mi \, \b r \. These 
"disk quantities" have been averaged using the "dcden" averaging 
(i.e. density-weighted average). Power-law fits for cases "f", "fdc", 
"dc", "ff 1 ", "eq", and "jet" are also considered but not shown (see 
section |3T4| for definitions of these other cases and more details). 

These "disk quantities" are fit for radii r dcden to r dcden , where 
r <fcden = X2 is typically set and /- dcdcn corresponds to N = 3 inflow 
times for time T? computed via Eq. 1 28 l except limited to r dcden < 



30r g to get an accurate fit across all models. A single inflow for time 
Tf is given as r dcden 



6.8 



. The stagnation radius (r s ), used in section 
is used to restrict r dcde " < r s . For some models (e.g. MB09D), if 
r s < I2r g or r dcden < 12r g , then we force r dcden = 12 and r dcden = 16. 

Table 1 10| also shows quantities M in — M H , M mw , and M w . For 
these wind quantities, fits are from r = 20r g up to = r dcdc " de- 
fined before, but > 40r g and r* < 100r g are enforced. 



For thick poloidal field models, roughly: p oc r 

\v,\ oc r- 8 , |w rot | oc r 03 , \b r \, Ifyl, \b\ oc r" 



- 6 , oc r-°\ 
, and M m , M„ oc r 13 . 
However, our A0.94BfN40 and A0.94BpN100 models have \v mt \ oc 
°. This occurs because our fiducial model (A0.94BfN40) was 
run for a longer time than any other models, while A0.94BpN100 
has no polarity inversions so also has accumulated magnetic flux 
over a longer time. So these models are choked out to a larger radial 
range. The flows are sub-Keplerian, with v Iot ~ 0.2v K at r ~ 10r g 
for the fiducial model (similar to in |Pen et al.|2003| >. 

For thick toroidal field models, roughly: po oc r~ 06 , p g oc r~ - 8 , 
\v r \ oc r~ 01 ', |w rot | oc r -03 , \b,\ oc and 1^1,^1 oc r _1 , M in ,M w oc 
r 13 _ r 2 w ith secular spin dependence. The flow is slightly sub- 
Keplerian in magnitude with v mt ~ 0.5uk at r ~ 10r g . 

For our thinner (TNM1 1) models with poloidal magnetic flux, 
roughly: p oc r~ 01 (depending upon initial /?), /? ? oc 9 , \v r \ oc 
(but depends upon spin), \v Iot \ oc r~ 03 , \b r \ oc f~ , |i>^|, \b\ oc r~\ and 
M- m ,M w oc r x -f- with secular spin dependence. These flows are 
mildly sub-Keplerian with v Iot ~ 0.6v K at r ~ 10r g for a/M = 0.99. 

The MB09 models are close to Keplerian both in profile (v TOt oc 
r -0 5 ) and value, as expected for weakly magnetized thinner disks. 

Consider an example usage of these fits. Table [4] gives 
Mj+ mw (r = 50r g ) ~ Mj + M mwo , and M w dominates. Table|lO|gives 
the radial power-law index, so the total mass ejected to large radii is 
Mj +mw ~ Mj+ mvl (r = 50r g )(r/50r g )" where n is from the second to 
last column in Table |10| For thick toroidal field models, the larger 
errors in M w (r) is because it flattens-out at larger radii. Perhaps 
steady-state has not been fully achieved by r ~ 100r g , or perhaps 
the self-similar region actually has a much shallower power-law 
index (consistent with M w oc r' 



Pangetal. 201 li 



In summary, our thicker disk MHD simulations lead to: flatter 
density profiles than the ADAF model, quite sub-Keplerian mo- 
tion, and flat pressure profiles at large radii due to the hot material 
reservoir. Our thinner disk models have flatter density profiles than 
the ADAF model and somewhat sub-Keplerian motion. Our thick 
poloidal/toroidal power-law results are stable to the chosen time- 
averaging interval, while some TNM11 models evolve to slightly 
steeper density profiles. 
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Table 10. Inner and Outer Radii for Least-Square Fits, Disk+Corona Stagnation Radius, and Fitted Power-Law Indices for Disk and Wind Flows 
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6.10 Resolution 

In this section, we determine whether our models are numerically 
converged. We consider both explicit convergence testing and so- 
called "convergence quality factors" that (e.g.) measure how many 
grid cells resolve various critical length scales (Sano et al. 2004 
|Hawley et al.|[20TT] [Shlokawa et al.||2012) |Sorathia et al.||2011| > 
Note that MRI quality factors like <2b,mri> 2<j,mri are not too useful 
for the new poloidal field models for which the MRI is suppressed. 

First, consider our overall resolution, box size, and numerical 
method choices. We knew that \m\ = 1 modes would dominate and 
even be required for (non-axisymmetric) accretion once the mag- 
netic barrier formed in poloidal field models. So, we chose to treat 
the poloidal and toroidal dimensions equally by ensuring the grid 
has a uniform grid cell aspect ratio, which also allows us to use less 
poloidal resolution to accurately capture the MRI lHawle y" et al.| 
|201 1| >. Also, because in = 1 modes generally dominate, Acp = 2n 
is required, as also found by Shiokawa et al. (2 012} and |Henisey| 
et al. ( 2009 1. In addition, other HARM-based GRMHD codes use 
the diffusive Toth scheme, while our code uses a staggered field 
scheme that treats field quantities more accurately. 

One indicator that an MHD simulation is unresolved is that 



MHD turbulence decays, magnetic field strengths drop, and fluxes 
secularly tend towards Novikov-Thorne values. These were used 
by |Shiokawa et al.|pOT2) and |Noble et al.| l |2010) as evidence that 
certain resolutions were insufficient to resolve the MRI. Except for 
the MB09D model, we see no evidence for such behavior. 

Just because turbulence does not decay does not mean the 
saturated state is converged. Convergence can only be ensured 
by performing explicit convergence testing and by resolving crit- 
ical length scales. a- mag > 0.4, <2„/„,, cor > 6, 2a.mri ^ 10, and 
2<».mri ^ 20 may be required to achieve tens of percent level con- 
vergence for the MRI (Ha wley et al.|20lT]|Shiokawa et al."| 2012l. 
These quantities are given in Table [8] as described in section [677] 

Consider how changing Nf changes the results for models 
A0.94BfN?c? and A0.94BfN30(r). Although initial p varies in 
those models, we already discussed how the magnetic flux near the 
BH has saturated independently from the initial magnetic flux. All 
tabulated results for these 3D models are similar to within ~ 30% 
fractional differences. Also, different realizations (A0.94BfN30(r)) 
and different fi (A0.94BfN100cl) are within ~ 10% fractional 
differences in Z7 H - Statistically, our averages are only accurate to 
~ 20% fractional differences, so values within this range are con- 
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sidered similar. However, we already saw from section |5T6| and sec- 
tion |5J] that there is significant power at higher |m|. While these 
higher \m\ modes only moderately affect the 6 - <p integrated quan- 
tities on the horizon or other radii, they significantly affect the time 
dependence of the solution. A plot (not shown) of (e.g.) efficiency 
vs. time shows more violent oscillations at lower N$. This suggests 
that the temporal behavior is qualitatively affected by how well- 
resolved the \m\ modes are, and only the N$ = 128 model shows 
temporal variability similar to the fiducial model. 

Consider how changing N r ,N e ,N,p by a factor of two 
changes the results for models A0.94BfN40 (fiducial model) vs. 
A0.94BfN100cl and A0.94BfN30(r). Changes in measured quan- 
tities are order tens of percent as discussed before. In addition, con- 
sider the azimuthal correlation length's m mode (m col , via Eq. {37}) 
for quantities p , u g , b 2 , u', b h 5,, F M , F™, A (and their absolute 
value versions) both in the "Disk" and "Jet". For A0.94BfN40, 
across all quantities, m COI ~ 6-14, except in the "Disk" we found 
w cor ~ 20 for b r and m cor ~ 15 for b 2 at r/r g = 4,8,30. On the 
horizon itself, where the disk is quite geometrically thin and the 
flow is causally disconnected from the rest of the solution, mag- 
netic field components in the disk+corona have m csx < 45. Even 
for field components, <2,„_ cor > 14 (via Eq. j44}, grid cells per cor- 
relation length) outside the horizon and Q mxm > 6 on the horizon. 
Beyond the horizon, the jet is always even better resolved than the 
disk. So, our lower-resolution choice of Afy = 128 is sufficient to 
resolve most structures beyond the horizon. 

Consider the vertical and radial correlation lengths in the 
"Disk". At r = ru,4r g , 8r s , 30r s , respectively, we find / P0 , C or ~ 
95, 58, 58, 57 giving 2/, C or,p ~ 3' 6> 6, 5 grid cells per vertical cor- 
relation length. Also, 42 c01 . = 108,37,27,17 giving <2/,corb- ~ 
3, 8, 12, 18. Also, the grid cells per radial correlation lengths are 
e».co,. P „ ~ 7,18,22,22 and Q ll£mb i « 7,16,21,22. Note that 
taking the spectrum of the averaged flow rather than the aver- 
aged spectrum leads to about twice higher apparent mode reso- 
lution for the vertical and radial correlation lengths. With ff 1 « 
0.06,0.13,0.29,0.59, this gives Ag, cos , p Jff' * 0.6,0.4,0.2,0.1 and 
^e.cor.b 2 /^ ~ 0-5, 0-7, 0.4, 0.3. Beyond the horizon, the jet is always 
even better resolved than the disk. Summarizing, this shows that the 
narrow density filaments are fairly resolved despite the strong mag- 
netic Rayleigh-Taylor instabilities, while the magnetic field that 
fills-in the region between the dense filaments is well-resolved be- 
yond the horizon. Across our poloidal field models, N s = 128 is 
optimal to resolve the compressed dense magnetic Rayleigh-Taylor 
filaments outside the horizon, while the magnetic field is marginally 
resolved even at N e = 64 (used for sweeping over spin). 

Consider the A0.94BtN10(HR, i.e. high-resolution) toroidal 
field models. The HR model gives a b and all 77's as quite similar. So, 
our lower-resolution toroidal field models are probably quantita- 
tively converged. Indeed, all our toroidal models have <2«.mri <: 10, 

2fl,weak,MRI ^ 10, <2(*,MRI ^ 20, <2<t,weak, MRI ^ 20, and ttmag ~ 0.4 

as required to well-resolve the MRI lHa wley et al.|[20TT) . The 
A0.94BIN10 model has ft, MR i > 10, £ e ,„eak,MRi > 10, Q<f, Mm > 58, 
G^.weakiURi > 32, and or mag ss 0.34. The A0.94BtN10HR model has 

Gfl.MRI 3 s 50, 2fl.wcak.MRI > 30, <2aMRI 3 s 170, <20, W eak,MRI > 80, 

and a mae « 0.38. (The stated Q's are limited by flow at r = r 
where 3 inflow times have passed.) So the MRI is probably well- 
converged. In addition, for A0.94BtN10HR, across all quantities 
(see list in previous paragraph) and locations (disk+corona and 
jet), the azimuthal correlation's m cor as 6-22 (typically ~ 10) at 
all radii r = rn, 4r g , 8r g , 30r g corresponding to <2,„, col > 12 (typi- 
cally ~ 20) grid cells per correlation length. Also, this corresponds 
to a typical azimuthal correlation length dif) COI ~ 0.9ff', so that the 



largest correlated azimuthal structures are about as extended as the 
half-vertical disk extent. These facts suggest that = 128 (all 
our toroidal field models have N$ 3> 128) is sufficient to well- 
resolve azimuthal structures. The vertical correlation lengths are 
resolved with g/.coi ~ 18-30 (typically ~ 25) cells across all quan- 
tities and all radii, while the radial correlation length is resolved 
with <2«,cor ~ 7, 20, 25, 25 grid cells at r = r H , 4r g , 8r g , 30r g , respec- 
tively, for both po and b 2 . In summary, our A0.94BtN10HR model 
is among the highest-resolved global MHD simulations. 

The 2D axisymmetric simulations are inappropriate for study- 
ing MCAFs (see also Igumenshchev 2009 ). Once magnetic flux has 
accumulated up to a saturation point even beyond the BH, accretion 
cannot occur in axisymmetry except through reconnection. Once so 
much magnetic flux has accumulated just beyond the BH, the en- 
tire flow rebounds backwards leading to low rjn. Eventually, mass 
builds up and forces magnetic flux back onto the BH leading to 
high 77h- Also, of course, 2D axisymmetric simulations cannot re- 
solve the non-axisymmetric MRI or sustain a magnetic dynamo. 

We now compare our resolutions with prior simulations of 
magnetic flux accumulation. Stehle & Spruit (2001) used a reso- 
lution up to N r x = 156 x 128 per decade in radius. So their 
simulations are roughly equally resolved to our fiducial models that 
have about half of the resolution per decade. They used van Leer 
interpolation, which has less than half the accuracy of our PPM- 
type interpolation. The 3D pseudo-Newtonian (PN) MHD simu- 
lations by Igumenshc hev et aL| ([2003 ) used a Cartesian grid with 
A(f> = n/2. Their inner-most cell size is 0.5r g at R m = 4r g . (a quite 
large _K in , see McKinney & Gammie 2002 ) Our fiducial model has 
dr ~ 0Ar g , dz ~ 0.037r g , and rsin0# ~ 0.097 r g at r = 4r g and 
has dr ~ 0.035r s , dz ~ 0.0076r g , and rsin8d(p ~ 0.033r g at r = m, 
so our z-resolution is about lOx higher. The 3D PNMHD energy- 
conserving PPM-type simulations by Igumenshchev (20081; Pun- 
|sly et al.|p009) are full A(p = 2n with N r xN e xN^ = 182x84x240 
for a comparable resolution per radii as our fiducial model, except 
very close to the BH where we have about 4x the ^-resolution. 



7 DISCUSSION 

Our simulations show that the accumulation of poloidal magnetic 
flux leads to a two-phase-like magnetospheric accretion flow that 
is dramatically different than the standard MRI-driven MHD tur- 
bulent accretion flow. The flow that develops in our simulations 
is conceptually similar to the "magnetically arrested disk" (MAD) 
flow ( TNarayan et al.|2003] >. While the standard weakly magnetized 
MRI-driven MHD turbulent flow has gas and magnetic pressures 
in force balance near the hole, the MAD state develops as mag- 
netic flux accumulates and magnetic forces balance the inflow's 
ram or gravitational forces. The originally-conceived MAD flow 
has a sharp magnetospheric boundary layer with a large density 
contrast at some radius, as confirmed by low-resolution 3D MHD 
simulations (e.g., fig. 13 in Igumensh chev et al.|2003| ; also seen 
in our 2D axisymmetric simulations). In these pioneering studies, 
accretion occurs primarily via diffusive reconnection events. 

Our high-resolution fully 3D simulations show that efficient 
non-axisymmetric magnetic Rayleigh-Taylor (RT) instabilities pre- 
vent the formation of the MAD's sharp magnetospheric barrier. 
Any additional magnetic flux that tries to accrete onto the BH is 
redistributed out in the disk by these instabilities. Also, we found 
that the magnetosphere geometrically compresses the dense inflow. 
We call this fully non-linear MAD flow a "magnetically choked ac- 
cretion flow" (MCAF), referring to the magnetic flux compressing 
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the dense inflow leading to enhancement of the magnetization over 
much of the horizon. Such a magnetic choke is analogous to chokes 
in man-made engines, within which it enriches the fuel mixture by 
partially shutting off the air intake. 

Our simulations confirm the brief 3D pseudo-Newtonian 
MHD simulations by Igumenshchev ( 2008 1 that also show MCAF 
formation. We also roughly confirm the magnetospheric QPO 
mechanism by |Li & N arayan (2004), which in their model drives 
some disk-based frequency at a vertical magnetic barrier. However, 
our 3D simulations show that the disk inflow interacts with the po- 
lar magnetic flux threading the rotating BH, which leads to a new 
"jet-disk" QPO (JD-QPO) mechanism based upon the BH rotation 
frequency. We also reaffirm that the BZ mechanism operates effi- 
ciently, except our low-spin thick-disk models do not form rela- 
tivistic jets due to mass infall. Otherwise, the BZ mechanism leads 
to powerful jets directly from the BH for our poloidal field models. 
We confirm the 3D GRMHD simulations by Tchekhovsko y et aL| 
( |2011| l, who showed that outflow efficiencies of 7/ > 100% are pos- 
sible once the BH with \a/M\ > 0.9 has reached poloidal magnetic 
flux saturation. As in other MHD simulations, the entire wind's out- 
flow rate is roughly M w oc r. 

We also confirm the results of Igumenshchev (2008 ) that 
a/M = MCAF models have low heat+outflow efficiencies of 
j] ~ few percent. One would expect heat+outflow efficiencies of 
77 ~ 100% even for a/M = if even the dense inflow were signif- 
icantly arrested (Nar ayan et al.|2 003). In particular, for a/M = 0, 
our thinner disk models have ~ 5%, while the NT efficiency is 
7/ NT ~ 6%. For the thick disk models, r) H < 0% and rj™ ~ 0%. In 
the simulations, the heavy disk inflow is relatively unmagnetized 
and not sufficiently slowed to achieve high 77. 

Radiatively efficient MCAF states (not studied in this paper) 
might still be hyper-NT efficient even for a/M = {Narayan et al.| 
|2003| >. Our thinner disk models have a high specific enthalpy such 
that n™ KE ~ 34% and 77™ ~ -28% for a/M = and ^ AKE ~ 64% 
and ~ -62% for a/M = 0.99. Emission of that free thermal 
energy would give a radiative efficiency of up to n ml i ~ 28% for 
a/M = and r] mi ~ 62% for a/M = 0.99. However, the trend 
with thickness is roughly t/ en oc 6' across our models, so we would 
predict that thin disks with < 0.03 (as relevant for BH x-ray 
binaries in the thermal state ; Kulkarni et al. 2011) would have 
no enhanced radiative efficiency. Thin radiatively efficient MCAFs 
should be studied with GRMHD simulations to check. 

We confirm the suggestion by Igumenshchev et al. (2003 1 that 
most prior MHD disk simulations used initial conditions that lim- 
ited the available magnetic flux. For example, a relatively radially- 
narrow torus can only have a relatively small amount of magnetic 
flux inserted if also keeping /? » 1 to allow for the MRI. Also, 
much of the matter and field can be ejected or remain beyond 
the torus pressure maximum rather than being accreted. After "Ph 
reaches a quasi-steady-state, one can test whether an MHD simula- 
tion is limited by such initial conditions. One computes 1 Ph/ v P ( , (i.e. 
[flux on hole] per unit [flux on the hole plus just outside the hole 
of the same polarity]) and also ^Vh/^s ([flux on hole] per unit [flux 
on the hole plus available within the stagnation radius for radial 
inflow]). One must compute both because there may be plenty of 
same-polarity magnetic flux beyond the hole, but it may not be ac- 
creting. One could also compute how much flux is available within 
the inflow equilibrium region. Both v P H / v P a ~ 1 and V P H / V P S ~ 1 for 
MB09D, so the initial conditions artificially limited the magnetic 
flux that can reach the hole. Also, it appears that 'Ph/'P, ~ 1 in the 
simulations by Beck with et al.|p009fr , who show much magnetic 
flux is ejected to (or remains at) large radii. 



A local ff-viscosity leads to a poor model of angular momen- 
tum transport for the simulations. The effective a is different than 
predicted by local stresses divided by either total or magnetic pres- 
sure for all our thick (H/R ~ 1) disk and poloidal field thinner 
(H/R ~ 0.3) disk models. The ff-disk theory only works well for 
our weak poloidal field models of either very thin (H/R ~ 0.05) 
or relatively thin (H/R ~ 0.3) disks. Large-scale magnetic confine- 
ment forces, rather than local stresses, may be acting on the dense 
inflow. Also, for our toroidal thick disk models, convection may be 
important because Reynolds stress dominates Maxwell stress. 

Interestingly, S rfM Ri ~ 1 (or S rfiM Ri ~ 0.25-0.5) gives the 
disk's saturated vertical field strength in our toroidal (or poloidal) 
field models. Also, the effective viscosity is a e g ~ 0.1-1 (a^n is at 
most 2x smaller for the new thinner models and at most 4x smaller 
for the new thick disk models), which is sufficiently consistent with 
observations (King et al. 2007]). 

MCAFs might explain observations of apparently highly effi- 
cient jets (fLaing et al.|1983]|Willott et al.|1999||Birzan et al.|20"04l 
|Ogle et al.|2006||Richards et al.|2006||Ghisellini et al.|2010||Mc>l 
Nama raet al.|201 l||Fer nandes et al.|2011||Punsly|2011[ >. However, 
more work is required to ensure the observations are accurately 
modelled. For example, the jet models in l Fernandes et al.| {201 \\ 
have factors (e.g. /, see |Willott et al.|1999| that can significantly 
change depending upon estimates of the invisible work done by 
the jet. Also, there are similar uncertainties in translating to the jet 
power in |McNamara et aL|l2011) as related to work done by bub- 
ble expansion (Birzan et al. 2004). Also, the Blazar estimates of jet 
power by Ghisel lmi et al.| (|2010l can be affected by models of the 
Doppler factor. 178Mhz observations by |Ogle eTa l. (2006) are af- 
fected by a lack of simultaneity between optical and radio emission, 
and the systems they observe could be low-luminosity radiatively 
inefficient systems that lowers their jet efficiency requirements. Our 
simulations may help constrain such jet power estimates. 

MCAFs may also help explain the variety of spectral states 
seen in BH x-ray binaries (Rutledge et al. 1999; de Gouv era"dal| 
|Pino & Laza rian 2005 ; Igumenshchev 2008 , 2009|. The accumula- 
tion of magnetic flux might lead to the low-hard, bright-hard, and 
hard:steep-power-law intermediate states with LFQPOs whose fre- 
quency is controlled by the "magnetospheric radius" (r,„) and with 
a persistent (mostly invisible) radio jet. The steep power law (very 
high) state could be due to reconnection as a field polarity rever- 
sal makes its way to the hole, where the reconnection with the jet 
on slightly larger scales makes the HFQPOs visible near the BH. 
Dissipation of the accumulated magnetic flux on the hole could 
cause the jet to light-up as a radio bright outgoing reconnecting 
plasmoid. The remaining disk without accumulated flux can lead 
back to the thermal (high-soft) state that has no jet and no (or very 
weak) QPOs. Notice that we distinguish between QPOs and jets 
being present vs. visible. Future work can test these speculations. 



8 CONCLUSIONS 

Black hole systems may have plenty of coherent magnetic flux 
available at large radii that can feed the accretion flow down to 
the black hole. Using fully 3D GRMHD simulations of extended 
radiatively inefficient accretion flows, we found that poloidal mag- 
netic flux readily accretes from large radii and builds-up to a natu- 
ral saturation point near the black hole independently from the ini- 
tial poloidal magnetic flux. The accumulated poloidal magnetic flux 
naturally leads to a highly non-axisymmetric "magnetically choked 
accretion flow" (MCAF), within which the MRI is suppressed. 



© 2012 RAS, MNRAS OOO.fTpH 



Magnetically Choked Accretion Flows 33 



In such a choked state, the polar magnetic flux forces accre- 
tion to occur through magnetic Rayleigh-Taylor instabilities and 
geometric compression. Near the black hole, the inflow is highly 
compressed, which leads to a highly magnetized state over most of 
the horizon. So the accumulated flux acts as an inflow choke (anal- 
ogous to air chokes in engineering that lead to fuel enrichment), 
which allows for greater than 100% jet efficiencies for \a/M\ > 0.9. 

Once the horizon's magnetic flux reaches a saturated state, 
the inflow pushes its way through the jet's bulging magnetosphere. 
This leads to a new jet-disk high-frequency quasi-periodic oscilla- 
tion (JD-HFQPO) mechanism driving spherical harmonic |m| mode 
oscillations with a frequency set by the black hole's field line an- 
gular rotation frequency at the jet-disk interface. So these HFQPOs 
may allow one to measure black hole spin. The coherence quality 
factor is highest in the jet and at the jet-disk interface (harboring 
large magnetic energy) rather than in the disk plane, so these HFQ- 
POs could be dominated by non-thermal emission. More work is 
required to test their observability. 

In contrast to our initially poloidally-dominated magnetic field 
models that develop the uncommonly-found MCAF state, our ini- 
tially toroidally-dominated magnetic field simulation results are 
similar to as seen in prior works. A significant new result is that 
such models do show large-scale dipolar magnetic flux formation 
near the horizon. However, only weak transient highly magnetized 
relativistic jets emerge. Higher resolutions or spins could promote 
this emergent dipolar field to launch persistent relativistic jets. 

As with any numerical study, more convergence testing is re- 
quired. Our toroidal field models satisfy the convergence quality 
factors of Qb.mri <: 10, Q^mri £ 20, and a mag a 0.4 in |Hawley etal.| 
l |2011) and r, 9, <f> correlation lengths are well-resolved. Our fiducial 
poloidal model well-resolves all correlation lengths. We explic- 
itly tested convergence with r, 9, (^-resolutions for a/M = 0.9375 
poloidal/toroidal models, and most values are converged to 30% 
and some values (e.g. T for poloidal models) are converged to 10%. 
We found that the convergence criteria of Hawley et al. (2011) are 
not generally sufficient (e.g. jet efficiency continues to change sig- 
nificantly in toroidal models) or applicable (e.g. MRI can be sup- 
pressed). Our code conserves energy during large-scale field recon- 
nection (as in |Lyutikov & M cKinney 201 1), occurring in magnetic 
Rayleigh-Taylor disruptions and field polarity inversions, but the 
reconnection rate may be set by grid diffusivity. 

In summary, the MCAF state and the spontaneous large-scale 
dipolar field generation seen in toroidal field models should be 
accounted for when seeking to understand SgrA*, M87, Blazars, 
high-powered quasars, black hole x-ray binaries, and other sys- 
tems. In future work, we will consider the spin dependence of the 
jet/wind power ( Tchekhovskoy & McKinney 2012) and of the JD- 
QPO, characterize the "magnetospheric radius" and determine its 
spin dependence, perform radiative transfer to identify the observa- 
tional signatures of MCAFs, and study the source of angular mo- 
mentum transport (e.g. turbulence vs. large-scale stresses). 
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Additional Supporting Information may be found in the online 
version of this article: Movie file. Movie of the fiducial model 
A0.94BfN40 showing the animated version of Fig. [4] (see caption 
alongside movie file for more detail). 



APPENDIX A: NUMERICAL METHODS 

Our simulations use the GRMHD code HARM based upon a con- 
servative shock-capturing Godunov scheme with 2nd order Runge- 
Kutta time-stepping, Courant factor 0.8, LAXF fluxes, simplified 
wave speeds, PPM-type interpolation (with no contact steepener, 
but with shock flattener based upon rest-mass flux density and total 
pressure) for primitive quantities (jpo, u g ,ff,ff\), a staggered mag- 
netic field representation, and any regular grid warping ( |Gammie] 
|et al.|2003[[Noble et al.|2006l|Tchekhovskoy et al.|2007) . 

One key feature of our HARM is its handling of grid cells 
with large magnetic energy per unit rest-mass energy. The code first 
tries to convert conserved to primitive quantities using the ideal 
MHD equations (Mignone & McKinney 2007). If that fails, the 
entropy conservation equations are tried. If that fails, the cold ideal 
MHD equations (with the failed grid cell having its internal energy 
averaged over its immediate neighbors) are tried. If that fails, the 
primitives are averaged over non-failed neighbors. The inversions 
are sought to machine accuracy. Density floors are used to avoid 
too low densities that lead to inversion failures. u g is chosen as to 
enforce « g /po < 50, then p is chosen as to enforce b 2 /p < 50 and 
b 2 /u g < 10 , and then these implied comoving changes are added to 
the stress-energy tensor in the ZAMO frame. Using a fixed (rather 
than comoving) frame avoids run-away energy-momentum gains. 

These inversion reductions and density floors mostly occur in 
the cold highly-magnetized jet with ordered magnetic field along 
which the flow moves. So, we set p = u g = in calculations 
(except in color images) if, e.g. for the "thick disk" simulations, 
b 2 /po > 30 to b 2 /po > 10 from rn to r = 9r g , respectively, as 
interpolated linearly with radius (tracing injected mass-energy, as 
in [Tchekhovskoy et al. 201JJ gives a bit more accuracy). Floors 
affect 5*, Q F only near the axes (not in the disk or collimated jets). 

Another key feature of our HARM is how the 3D polar axis 
is treated ( |McKinney & B landford 2009). To resolve the stag- 
gered field component 5 (2) , a machine error polar cut-out of size 
d9 = 10~ 13 is used. The values of B} 2> on the polar cut-out are 
copied asymmetrically across the pole. All other primitives are in- 
terpolated from active cells to ghost cells at the same physical loca- 
tions and also to the polar cut-out. All fluxes are set to vanish at the 
pole, except for the induction update to ■\J-gB} 2} . The small 9 fac- 
tors in y/^g cancel on both sides of the induction equation allowing 
one to recover 6 (2) on the polar cut-out. This procedure introduces 
a complex phase to the otherwise singular polar axis, where for 
general discontinuous flows the only constraint is from V ■ B = 
that forces B {2) to be continuous along 9 across the pole for each <f>. 
This transmissive polar boundary condition allows the flow to pass 
through the polar region with only additional dissipation due to the 
polar triangular mesh as compared to a more uniform mesh. 

Another key feature is the connection coefficients are cor- 
rected to obtain zero body forces to machine error for constant 
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pressure. Nominally, constant pressure flux and source terms do 
not cancel, and this secular differencing error sits at the polar axes. 
For constant pressures, the relevant equation of motion is 

d,( V^sO = -dj( ^6[) + yFgfifc. (Al) 

Assuming d,{yp-g) ~ and dj(p tot ) ~ 0, then = -d v (-\f-g) + 
■\f—gT* k or T k VK = d v (-\f zr g)/ yf—g. We need to subtract off face- 
related flux-positioned values and add center-related conservative- 
positioned values. Let a grid cell indexed by i,j,k have grid face 
position xt = x!i{i, j, k) and upper grid face position of 4 = 
Xj(i + dx"6 [ K ,j + dx K S 2 K ,k + dx K 6 > K ) (k chosen, not summed) such 
that the cell size is Ax* = x\ — jci and cell center position is 
j<£ = {x% + xV)l% First, compute C K = Fj^ as using a continuous 
approximation to the derivatives at x c . These are computed using a 
modified dfridrQ in Press et al. ( 1986) to achieve a uniformly low 
error. The conversion from internal to physical coordinates is gen- 
erated from a simplified continuous approximation to the derivative 
with a difference size of \.0~ 5 dx K , which depends upon range in x" 
as necessary to ensure the differencing scales with the grid. Second, 
compute D K = r*^ using discrete differences across the actual grid 
cell via D K = ( V^?U+] - -\FI[ x -])/( Ax * ^FiUcY)- Now, chang- 
ing each k term independently can lead to too large changes in some 
limits. Instead, we construct the weights S K = 1CT 300 + |r^| (sum 
over n) and let dS K = D„- C K , W m = ir^l/S* (no sum over k or fi). 
Finally, for each ji and each k, the correction is 

AT* = dS K W K , (A2) 

(no sum over k or fi). The new connection gives machine error can- 
cellation between source and flux differencing of pressure. 
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